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Column Behavior in the Plastic Stress Range 


JOHN E. DUBERG* ann THOMAS W. WILDER, IIIt 
Langley Aeronautical Laboratory, N.A.C.A. 


SUMMARY 


This paper presents a summary of the significant findings of a 
theoretical study of column behavior in the plastic stress range 
When the behavior of a straight column is regarded as the limit 
ing behavior of an imperfect column as the initial out-of-straight 
ness approaches zero, the departure from the straight configura 
tion occurs at the Shanley tangent-modulus load. Without such 
aconcept of the behavior of a straight column, one is led to the 
unrealistic conclusion that lateral deflection of the column can 
begin at any load between the tangent-modulus value and the 
Euler load, based on the original elastic modulus. 

The behavior of a column with vanishing initial out-of-straight 
ness at loads beyond the tangent-modulus load depends upon the 
stress-strain curve for the material. A family of load-lateral 
deflection curves is presented for idealized H-section columns of 
various lengths and of various materials that have a systematic 
variation of their stress-strain curves. These curves show that 
for columns whose material stress-strain curves depart gradually 
from the initial elastic slope, which is characteristic of stainless 
Steels, the maximum column loads may be significantly above 
the tangent-modulus buckling load. If the departure from the 
elastic curve is more abrupt, such as for the high-strength alumi 
hum or magnesium alloys, the maximum load is only slightly 
above the tangent-modulus load 


INTRODUCTION 


i RECENTLY, it had generally been accepted 
that the correct theory of column failure in the 
inelastic range of stress was the double-modulus theory. 
This theory predicts that the load at which bending 
starts and the maximum load that a column can sup- 
port are the same and can be obtained from the Euler 
equation 


P = WEI/L* 


by substitution of a reduced modulus for Young's 
modulus. The reduced modulus is obtained by as- 
suming that, at the start of bending of an originally 
Straight column, the direction of straining of the ele- 
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ments on the convex side of the column reverses. When 
such reversal of strain occurs in the plastic range, in- 
crements of stress are related to increments of strain 
by the original elastic modulus. The amount of the 
cross section over which the strains reverse is determined 
by the condition that there shall be no change in load 
during the bending process. In 1947, Shanley' was 
able to show for a simplified column that, if the load is 
allowed to increase during bending, bending of the 
column can start at a lower load than the reduced- 
modulus load. The load for which he showed this to be 
true was the tangent-modulus load and can be obtained 
from the Euler equation by substitution of the tangent 
modulus in place of Young’s modulus. Shanley drew 
conclusions concerning the behavior of columns on the 
basis of the behavior of the simple model and certain 
experimental observations. 

To clarify the behavior of columns in the plastic 
range, a study was made with the following threefold 
purpose: 

(1) To establish a definition of the critical load for 
a real column. 

(2) To study the mechanism of column action be- 
yond the critical load. 

(3) To establish the relation between the maximum 
column load and the stress-strain curve for the material. 

In order to make this study, two models were chosen. 
that is, one 
second, an 


Cne was a spring-supported rigid column 
that has a concentrated flexibility—the 
idealized H-section column that has its flexibility dis- 
tributed along its length and consists of two concen- 
trated flanges separated by a web of negligible area but 
of infinite shear rigidity. This study is reported in 
sone detail elsewhere,” and only the more important 


results are summarized here. 


SYMBOLS 


= width of column 
lateral deflection of top of spring-supported column 


(dy excluded) 
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dy = initial lateral deflection of top of spring-supported 
column 
e€ = strain, unit strain, or displacement 
e = critical spring displacement 
E = Young’s modulus 
F = force 
I = moment of inertia 
ky = stiffness of fixed-end spring 
ke = stiffness of ratchet-ended spring 
E = measure of length 
n = exponent used in analytical form of stress-strain curve 
P = load on column 
P, = load on straight column causing a stress of o; or a dis 


placement of e; 
p = radius of gyration 
= mid-height deflection of H-section column 


= stress 
v1 = stress appearing in analytical form of stress-strain 
curve 
Subscripts 
e = Euler load 
max. = maximum load 
RM = reduced-modulus load 
r = tangent-modulus load 


THE CRITICAL LOAD 


In the elastic range of stress the critical load for a 
column—that is, the load at which bending starts 
is unique and is given by the well-known Euler formula. 
The situation is not so simple in the inelastic range of 
stress, however, and the source of the difficulty lies in 
the character of the stress-strain relations in the plastic 
range. In the inelastic range of stress, at least for uni- 
axial states of stress, increments in stress are related to 
increments in strain by the tangent modulus of the 
material, but decreases in stress are related to strain 
by the original elastic modulus. 

The flexural stiffness of a column in the bent posi- 
tion, upon which the strength of the column depends, 
is a function not only of the load but also the amount 
of strain reversal which has taken place. Thisleads toa 
lack of uniqueness in defining a critical load. 

In order to demonstrate more clearly this lack of 
uniqueness in defining a load at which bending starts, 
the simple spring-supported column shown in Fig. la 
was analyzed. The inverted tee is rigid and free to 
rotate and translate vertically about the intersection of 
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Fic. 2. Column behavior at critical loads. Fig. 2a (lef 
Initial slopes of load-deflection curves. Fig. 2b (right). _Initig) 
stress and strain changes at center cross section. 


the cross of the tee. The column is supported at each 
end of the horizontal leg by identical sets of two elastic 
springs. One spring is fixed at the far end. The other 
spring has a ratchet attached to its far end which per- 
mits no additional strain in the spring as the displace. 
ment of the end of the horizontal member exceeds ¢, 
If the end of the horizontal member moves upward 
after exceeding e,, the ratchet immediately catches. 
The combined load-displacement relation of each spring 
system is as shown in Fig. lb. This relation may be 
regarded as a simple stress-strain curve that includes a 
plastic region and the phenomenon of strain reversal. 
If such a column is assumed perfectly straight and its 
load is free to change during bending, an infinity ol 
loads can be found at which the top of this column can 
start to assume a deflected position. This range o/ 
loads is included between the tangent-modulus load for 
this column and the Euler load. Those loads between 
the tangent-modulus load and the reduced-modulus 
load require an increase in load during the bending 
process, whereas those loads between the reduced 
modulus load and the Euler load require that the load 
decrease during the bending process. 

The actions of the spring system at the start of bend 
ing can be described by locating the instantaneous 
center of rotation of the rigid column. At the tangent 
modulus load the center is at the end of the horizontal 
leg on the side that is not loading, so that no dis 
placement of that spring system occurs. For loads 
between the tangent-modulus load and the Euler load 
the center lies on the horizontal leg between the two 
ends. At the Euler load it is at the opposite side, and 
no displacement of that spring system occurs. It 's 
evident then that for all these loads reversal of the 
spring system on one side always occurs except at the 
tangent-modulus load when there is no change in the 
displacement of that side. 

An analysis made for the idealized H-section reveals 4 
similar but somewhat more complex behavior of a per 
fectly straight column (see Fig. 2). There is the same 
range of loads at which a perfectly straight plastic col: 
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mn can start to assume a deflected shape. Those 


u 
loads between the tangent-modulus load and the re- 
d-modulus load are associated with an increase in 


Those between the reduced- 


duce 
load during deflection. 
modulus load and the Euler load are associated with a 
decrease in load. The essential difference in the col- 
ymn actions for this range of loads is the pattern of 
strain reversal which occurs during bending. At the 
tangent-modulus load, none of the strains over the 
whole volume of the column reverse at the start of 
The strain, however, is just stationary at 
At the 
reduced-modulus load, the reversal is complete over one 
side of the column, just as it is usually presented in the 


bending. 
the center of the column on the convex side. 


double-modulus theory, and the strains are stationary 
along the line that separates the reversed and un- 
reversed regions. At the Euler load, all the strains re- 
verse over the entire vglume of the column except for 
one point at the center of the column on the concave 
side, and there the strain is stationary. Because at 
these three loads the distribution of stiffness over the 
length of the column is the same, the instantaneous de- 
flected shape is a single sine wave. 

Knowing the deflected shape of the column at tan- 
gent-modulus and Euler loads and realizing that the 
instantaneous center for strain of the center cross sec- 
tion of the column is on the convex side at the tangent- 
modulus load and on the compression side at the Euler 
load, it is possible to compute the initial slope of the 
load-center deflection curve for the columns. The re- 
sults obtained here for the H-section are equally valid 
for all columns of constant symmetrical cross section. 

At the tangent-modulus load 


dP/dd = (b/2p*)P7 
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Fic. 3. Four force-displacement combinations for the spring 


systems. 
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Load-deflection curves for spring-supported col 
Fig. 4b (bottom). Detail of Fig. 4a. 


Fic. 4a (top). 
umn with initial deflection. 


and at the Euler load 


dP/dd = —(b/2p*)P. 


The quantity 6 is the distance between the extreme 
fibers in bending, and p is the radius of gyration of the 
cross section. At the reduced-modulus load, this slope 
is zero. 

The analysis of a perfectly straight plastic column 
leads to a range of critical loads, but only one can be 
No real 
column is perfectly straight, and therefore it is more 


accepted as having meaning for real columns. 


reasonable to define the critical load as one based on the 
behavior of a slightly bent column as the initial lack of 
straightness vanishes. In order to demonstrate the re 
sults of such a point of view for defining the critical load, 
a more complete analysis of the spring-supported col- 
umn was made in which a small eccentricity, do, of the 
top of the column was included. In making such an 
analysis, it is necessary to keep track of the displace- 
ments of the individual spring systems to be certain that 
the correct force-displacement relations are being used. 
The possible combinations for the spring system are 
shown in Fig. 3. 

The results of the analysis of the spring column with 
various amounts of initial eccentricity are shown in 
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Also shown are the limits of the regions in which 
The most 


Fig. 4. 
the various force-deflection relations exist. 
significant point to be observed is that, as the initial 
imperfection decreases, the tangent-modulus load is 
the limiting load at which there is a sustained increase 
in the bending deflection of the column. It should also 
be noted that the reversal of the spring deflection al- 
ways occurs below the tangent-modulus load and occurs 
just at that load as the initial imperfection vanishes. 


THE MaAximuM Loap 


The load-deflection curves for the simple column, 
given in Fig. 4a, are all approaching the reduced-modu- 
lus load at large deflection. This is a consequence of 
the linearity of the spring systems and is to be expected 
if one considers the deflected column with reversal as a 
new elastic column whose stiffness is measured by 
the reduced modulus. This is not typical of the ac- 
tions of real columns, because as the column loads, 
the strains increase and there is a continual reduction of 
the tangent modulus. How rapidly the tangent 
modulus decreases depends on the shape of the stress- 


strain curve. Therefore, it is to be expected that the 
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maximum strength of a column of a given geometry wij) 
depend on the stress-strain curve of its material. 


To study the effect of the shape of the stress-strain. 


curve on the maximum strength of a column, an ana. 
lysis was made of the behavior of the H-section ¢oj. 
umn after it had become critical and started to bend at 
the tangent-modulus load. It was assumed that the 
analytical stress-strain curve for the material of this 
column was of the form suggested by Ramberg and 
Osgood.* This form is summarized in Fig. 5. The 
stress a; is usually close to the yield stress of the mate. 
The most significant parameter for this study js 
Low values of 1 correspond to gradu. 


rial. 
the exponent 7. 
ally curving stress-strain curves, and, as increases, 
the curvature changes more rapidly at the knee, 
Values of 1 in the neighborhood of ten are associated 
with aluminum alloys, whereas values about half this 
apply to stainless steels. Magnesium and the low-car. 
bon steels have stress-strain curves that correspond to 
values of of 30 or greater. 

The study of strain history of the flanges of a column 
that starts to bend at the tangent-modulus load shows 
that the same general strain history occurs regardless 
of the shape of the stress-strain curve. This is shown 
pictorially in Fig.6. At the start of bending, reversal o 
the stresses begins at the center of the convex side of th 
column. As the deflection increases, the region 0 
reversal spreads rapidly over the convex side of the 
column and, at maximum load, is complete over the 
whole convex flange. After the maximum load is 
passed, reversal of strain spreads into the concave side 
of the column. 

A summary of the load-deflection results that were 
obtained for the H-section column when the stress 
strain curve and the tangent-modulus load were system 
atically varied is given in Fig. 7. The loads are given 
in terms of a load P; that produces an average stress ¢ 
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the straight column. Two results are significant: 
First, the smaller the value of m, the higher the maxi- 
mum load in relation to the tangent-modulus load; 
cond, the smaller the value of m, the greater the de- 
fection at which the maximum load occurs. The re- 
sults obtained for the maximum load are summarized 
intwoforms. In Fig. 8, the ratio of the difference be- 
tween the maximum load and the tangent-modulus load 
to the difference between the reduced-modulus load 
(had the column remained straight) and the tangent- 
modulus load are plotted as a function of tangent- 
modulus load. The maximum load of columns critical 
in the plastic stress range exceeds the tangent-modulus 
load by a fairly constant percentage of the difference 
between Pry and Pr. In Fig. 9, the ratios of the 
P... to Pr are plotted as functions of the tangent- 
modulus load. In the plastic range the percentage 
increases in load over the critical load are roughly con- 
stant and are greater the smaller the value of n. 

The dashed portions of the curve in Figs. 8 and 9 have 
no practical significance. Their shape is a consequence 
of the fact that the Ramberg-Osgood stress-strain curves 
have no proportional limit and are therefore not linear 
in what would normally be the elastic range. A more 
orrect interpretation would be to consider the ordinates 


obe unity in the dashed regions. 


CONCLUSIONS 


(1) If the behavior of a perfectly straight column is 
regarded as the limiting behavior of a bent column as its 
initial imperfection vanishes, then the tangent-modulus 
load is the critical load of the column—that is, the load 
at which bending starts. 

(2) It is possible to compute the maximum load a 
column can support, but in order to do so it is necessary 
to include a proper treatment of the phenomenon of 
strain reversal. 

(5) In the plastic range the maximum load a col- 
umn can support is larger than the tangent-modulus 
load and, for a given column, is greater the more 
gradual the change in the slope of the stress-strain 


curve in the neighborhood of the yield stress. Prac- 
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tically, this means that columns made of stainless steel 
can have maximum loads significantly above the tan- 
gent-modulus load, whereas columns made of the high- 
strength aluminum and magnesium alloys cannot. 
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Second Approximation to Supersonic 
Conical Flows’ 


FRANKLIN K. MOORE? 
Cornell University 


SUMMARY 


The method of expansion of the velocity potential in powers 
of a small geometrical parameter is applied with a view to find 
ing a second-order correction to the linear theory for steady, 
isentropic, supersonic flow of an inviscid perfect fluid about a 
body having conical symmetry. A body lying entirely within the 
tip Mach cone is treated in detail. 

A method is developed for formulating, at a mean plane, bound 
ary conditions at a thin body, and the particular integrals of the 
differential equations applicable to the second-order velocity per 
turhations are developed in integral form. The technique given 
by Lighthill is applied to determine the upstream boundary con 
ditions consistent with the presence of an attached conical shock 
wave. The second approximation to the pressure coefficient on 
an “arrow-head”’ airfoil at zero angles of attack and yaw and 
which is infinite in downstream extent is computed by this tech 


nique. 


INTRODUCTION 


— LINEARIZED THEORY of supersonic flow in an 
inviscid perfect fluid has reached a high state of 
development. Because of its approximate nature, 
however, this theory is subject to error when applied to 
airfoils of finite thickness or angle of attack. 


no simple way of estimating this error, except in those 


There is 


few cases for which higher approximations or exact solu- 
tions are known. The results of the linear theory for 
the pitching moment of airfoils are particularly sus- 
pected. 

It would be valuable, therefore, to improve the linear 
theory for three-dimensional cases. This problem is 
approached in this paper via the method of expansion 
of the velocity potential in powers of a small geometri- 
cal parameter (thickness, say). This method will be 
applied to steady, isentropic problems of thin airfoils 
having conical symmetry—i.e., cases wherein the free 
stream is disturbed by a body generated by rays from a 
common focal point. 

If a velocity potential is expanded in powers of a 
thickness parameter and the resulting series is intro- 
duced into the equations of isentropic compressible 
fluid motion, a linear hyperbolic differential equation 
similar to the wave equation is obtained from terms of 


first order in thickness for supersonic velocities. The 
Received August 5, 1949. Revised and received April, 1950 
* Based on a Ph.D. thesis presented to the Graduate School 
of Cornell University in June, 1949. The author acknowledges 
the guidance and assistance of Profs. W. R. Sears and Y. H. Kuo, 
of the Graduate School of Aeronautical Engineering, Cornell 
University. 
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solution of this equation, satisfying the boundary 
conditions on the obstacle, formulated to the first order 
is called the first approximation. If terms of second 
order are collected, a linear differential equation results 
which involves the first approximation and whose soly 
tion, subject to the proper boundary conditions, repre 
sents the second approximation. 

Certain aspects of this problem require discussion 
prior to any detailed analysis: 

(a) Isentropy and irrotationality: The problems to 
be considered will be ones in which dissipative phe 
nomena may appear only within the bow wave (and 
Since the 
discontinuous entropy change across such a bow shock 


then only if the bow wave is a shock wave). 


wave is of third order in thickness, it may be neglected 
in the development of a second-order theory. Hence, 
the entire flow field may be considered isentropic and, 
consequently, irrotational. One is then justified in 
introducing the velocity potential. 

(b) The conical nature of the first and second approxi 
mations: It has been shown by Lagerstrom! that any 
exact flow over a body having conical symmetry will be 
conical (i.e., flow properties such as velocities and pres 
sures will be constant along any ray from the focal 
point of the body), provided only that viscosity does 
not act. This, of course, assumes the existence of the 
solution. A nontrivial solution does not exist for sub 
sonic flow over an infinitely extended conical body. 

The solution exists for supersonic flows, due to the 
“rule of forbidden signals.’’ Since disturbances maj 
not propagate upstream, a body, such as an “‘arrow 
head”’ pointed into the wind, may be terminated to pro 
vide a physically possible flow without affecting the flow 
upstream of the termination. The resulting flow may 
then be considered the consequence of infinitely e 
tended boundary conditions with no sacrifice of physical 
reality. 

Thus, one concludes that, since the exact solution 1s 
conical and each term of the proposed expansion of @ 
velocity, say, is of different order of magnitude, each 
term of the expansion must be conical. 

(c) Justification of the expansion of the velocity po 
tential in powers of a thickness parameter: One expects 
a solution to the differential equations of motion to be 
analytic in terms of the thickness parameter near the 
surface of a thin airfoil. If this is correct, the solutions 
can be expanded in power series in this parametel. 
Of course, local breakdowns are expected in the neigh 
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porhoods of stagnation points, discontinuities in body 
curvature, or singularities of the differential equations. 

The clearest indication of the correctness of this pro- 
cedure is provided by consideration of the boundary 
conditions at the body. For a flat airfoil, the first- 
order solution must be of order ¢, where ¢ is the thick 
This is dictated by the boundary 
On the air- 


ness parameter. 
condition on normal velocity at the airfoil. 
joil, except at or near singular points of the first-order 
solution (e.g., leading edges), the same boundary con- 
dition then requires that the next nontrivial approxi- 
mation be of order e*, and soon. This statement may 
be checked against the boundary condition formulated 
in Eq. (29). Thus, the power expansion in ¢ is the only 
expansion which will satisfy the boundary conditions on 
the body. Recent work by Lighthill’ indicates that, 
to second order, the power series expansion in € prop- 
etly describes conditions at the Mach cone as well. 

Further, it is expected that the solution, to each 
order, will be analytic at some plane (e.g., a plane of 
symmetry) near the thin body; thus, the value of each 
approximation on the body can be expressed as a power 
series in e, the leading term being its value on the plane 
ofsymmetry,say. As will be shown later, one may then 
formulate boundary conditions at the plane of symi- 
metry to all orders. 

For purposes of this paper, the above statement need 
not be rigorously proved. One will need only to sup- 
pose that, if the first-order solution is taken to be of 
unit order, its value on the body differs from that on 
the plane of symmetry by a quantity of order «. This 
can be verified for any particular problem. 

Broderick* ‘ and Lighthill’ have solved correspond- 
ing problems of the supersonic flow about bodies of 
revolution and, in order to satisfy boundary conditions 
at the surface, have found it necessary to introduce 
nonregular terms involving the logarithm of the thick- 
ness parameter. This is a consequence of the fact 
that the axis of the body (corresponding to the plane 
of symmetry in the present case) is a singularity of their 
solutions. 

(d) The shock wave: In the case of supersonic conical 
flows over slender bodies, an attached conical shock 
wave may be presumed to occur in the vicinity of the 
tip Mach cone. In the linearized theory of such flows, 
itiscustomary to assume that the shock strength is zero 
and that the shock wave is located at the tip Mach 
Lighthill? has pointed out that the region near 
the shock wave requires special investigation in this 
He has worked out a technique for treating 
this region and has shown that (1) the assumption of 
zero shock strength is valid for the first (linearized) 
approximation, and (2) the true shock strength is of 
He presents 


cone. 


theory. 


second order in disturbance velocities. 
formulas that provide a quantitative second-order 
description of the flow near the shock wave in terms of 
results of the linearized solution. Lighthill's technique 
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will be used herein to provide boundary conditions at 
the Mach cone. 

(e) The pressure Since the 
namic forces on the body are of primary interest, the ob 
jective of the analysis to follow will be to calculate the 


coefficient; aerody- 


pressure coefficient on the body, which may be written 


as follows: 


C= ~= CO) 4+ 2CO 4+... (1) 


The velocity vector is 


ifU + eu? + eu? + 1+ jlev? + 


k[lew +... 


the 1 vector being in the streamwise direction, and 


B=V M2 — 1 (: 
whence, applying the isentropic equations of motion, 
U) (4) 


CY = —2(u 


U)? + (@®/U)? + 
w/U)?] (5) 


C® = —[2(u™/U) — Bu, 


Thus, one needs all three first-order cartesian velocity 
components and the second-order correction to the 
streamwise velocity component, only. 


NOTATION 


a 2 = cartesian space coordinates 
n = conical radial coordinate 
= Tschaplygin radial coordinate 

w = space, conical, and Tschaplygin angular co 
ordinate 

ae = complex coordinates in the Tschaplygin plane 

e = pressure 

p = density 

e = pressure coefficient 

U = free-stream velocity 

M = Mach Number 

B = a constant based on free-stream Mach Number 

u,v, % = cartesian velocity components in the directions of 
x, y, and z, respectively 

e, A = angles—-see Fig. 2 

¢ = velocity potential 

a = velocity of sound 

7 = ratio of specific heats 

Vv? = the Laplacian operator 

iL“ = geometrical quantities in the Tschaplygin plane 


see Fig. 3 
v,f,¢,h = flow functions 


A,D,K = 


Subscripts 


constants 


© signifies evaluation in the free stream 

» signifies evaluation on the body. 

O signifies evaluation on the plane of symmetry. 

p signifies particular integral. 

¢ signifies complementary function. 

The subscript notation for partial differentiation is used where 
convenient. 
Superscripts 

(n) denotes the order of approximation 

Primes denote ordinary differentiation. 
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Applying this transformation to Eqs. (9) and (10), it; . 
WwW TIP MACH CONE found that ( 
f Vu (s, 0) = Va = Vw) = (12 a 
Z — 2 Ss con 
Vu’? = 8 (1 <2)3 [2h, + sl — S*)itee) (13 in t 
solu 
vor = = : — [(1 + 4s? — s*) cos wh, + task 
(i — $7)" exal 
s(1 — s*) cos wh,; — (1 — s*)? sin wh,y] (14 
| 4 
Vow = — Fn [Cl + 48° — 54) sin wh, + Si 
s(1 — s4) sin why, + (1 — 5)? C08 why] (15) | Po 
lorn 
where the 
h = A®[_D*u®” + yD? 4 wn? (16 _— 
fied 
Fic. 1. Diagram showing notation and Mac 
¥-1 M.! wh 
THE GENERAL CASE D?=1+ _ M.?; A?=2 Ve (17 valu 
(1) Differential Equations for the Velocity Potential - mett 
cigar E ’ The following transformation may now be introduced Tl 
For the scheme shown in Fig. 1, the isentropic equa- : 
tions of motion are case’: —§=se* Is actus 
shocl 
(a2 — 9, )err + (a* — ¢y7)Oyw + (a? - ¢,2)¢,. — and one finds that Macl 
2629, Pry — 2yPHyz — 26292: = 0 (6) V2 = 4(d2/a¢d8) 19 = 
9 ) 9 9 9 = athe 
at = a.’ + [ly — 1)/2] (U* — on* — @%* — 924) WY) which makes it possible to define the following fune ordin 
It is assumed that tions: form 
g= Uz + og” + Ce” +... (8) AS) + gi(¢) DAu" 
folf) + gelO) Ay" (20 
Introducing Eqs. (7) and (8) into Eq. (6) and collect flo) + £3(0) Aw" 
ing terms of order ¢€ and e’, respectively, one finds ; oe where 
a a a (9) and to — a aap tor % ae wtepiaeaaass: mine 
ul 22 \ integrations by parts in Eq. (13), as follows: appre 
Brrr — gy? — ¢2 = wt «+ @ . ulatec 
= sa? = Deyo" gV* 4 gabe] — The : 
_ = - (3 + Z 9 u.*) gr? + gy P? + o:0*| R Zi-— #as | + + poh 
Eq. (9) is the equation of the linear-perturbation _— + + om ia ee | 
theory. / g, dé | if = (21 R<) 
The solution for ¢’ will be composed of a particular ee Lio! 
solution of Eq. (10) and a complementary solution Expressions for v,‘? and w, similar to the above can Pi 
satisfying the homogeneous counterpart of Eq. (10), be obtained from Eqs. (14) and (15). flow u 
which solutions, taken together, satisfy the boundary It is clear at this point that the solution given in Eg hacia 
conditions appropriate to the particular case considered. (21) will be infinite at the tip Mach cone (s = 1) unless at R 
: ; the derivative in the first right-hand term vanishes sileed 
(2) Teensformations and Particular Integrate there. This is an example of singular behavior that 
A transformation introduced by Tschaplygin® is ap- concerned Lighthill in the investigation already met 
plicable to conical flows and changes the three-dimen- tioned. In the cases considered in detail here, in which 
sional hyperbolic equations for velocity components the body lies entirely within the tip Mach cone, this 
into two-dimensional Laplace's equations. The trans- derivative does in fact vanish, since the linearize heoons 
formation is introduced as follows: theory leads to vanishing perturbation velocities oF contin 
- = this cone. Consequently, the second-order solution 's aad he 
7 yf. =e ts 2% ” free of this particular difficulty, and one needs Light | jy. 
n=B w=tan-!-; 5 I y» hill she 


x y 1+ V1l—y, _ hill’s technique only to establish the upstream boundary | Rat 
(11) condition. i 
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The Boundary Condition on the Body 


3 


One prescribes simply that the flow through the 
body surface vanish to the second order of approxima- 
tion. It would be convenient to satisfy this boundary 
idition at some mean plane, as is customarily done 
in the linear theory. The assumption of analyticity of 
solutions at this plane enables one to do this. The de- 
tails of this procedure will appear later in an illustrative 


col 


example. 
4) The ‘‘Boundary Condition” at the Mach Cone 

Since it is intended to deal herein with velocity com- 
ponents rather than the potential, it is necessary to 
formulate boundary conditions on these quantities at 
the Mach cone. As has already been mentioned, the 
formulation of correct boundary conditions to be satis- 
fied by the perturbation velocity components at the tip 
Mach cone is based upon Lighthill's detailed analysis’ 
of the conditions near the shock wave. It may be of 
value to review briefly the salient features of his 
method. 

The correct upstream ‘“‘boundary conditions’ are 
actually the shock-wave equations, applied at the 
shock-wave location, which is upstream of the tip 
Mach cone and is therefore outside the region of appli- 
cability of the method set forth above. To obviate this 
difficulty, Lighthill introduces a transformation of co- 
ordinates from 7, w to R, w, which may be written in the 
form 

n= Rt er(w) + €r(w) + 
where the successive functions 7‘"’(w) are to be deter- 
mined as the calculation progresses. The method of 
approximations in successive powers of ¢ is then form- 
ulated in terms of the new independent variables R, w. 
The singularity of the differential equations is thus 
shifted to the location R = 1, instead of 7 = 1. The 
choice of the functions r‘"(w) is made, however, in 
such a way as to avoid divergence of the method at 


R= 1. Itis then found that the shock wave occurs for 
R < 1--i.e., in the region of convergence. 
Lighthill next introduces the Rankine-Hugoniot 


shock-wave equations to connect the conical flow to the 
flow upstream of the shock. This leads him to fictitious 
boundary values assumed by the velocity components 
at R = 1. In particular, consider the perturbation 
velocity potential 
g = Uxf(R, w) 
= Uxflef(R, w) + &f(R, w) + | 


According to the shock-wave equations, this quantity is 
i.e., it has the value zero 
Light- 


continuous across the shock 
just behind the shock in the present problem. 
hill shows that the corresponding values to be applied at 
R = | are 


FLOWS 331 


fY1,w) = f™(1,0) =0 
Sr (1,0) = 0 
2 ¢ (22) 
1M . fr 
fy, w) = y+ lim - 
, 2 p> | k-1V1—-—R 


Furthermore, in any case where the body lies entirely 
within the tip Mach cone, it is found that r(@) is 
identically zero. This means that, at any point in the 
flow field, » and R will differ only to second order in e. 
Thus, away from the vicinity of the shock wave, the 
independent variable R in the functions f'(R, w) may 
be directly replaced by y without introducing any error 
greater than order e’*. 

It follows from the fact that r‘''(w) is zero that the 
differential equations satisfied by f"(R, w) and f(R, 
w) are just the equations equivalent to Eqs. (9) and 
(10) and, hence, to Eqs. (12)—(15) (the only differences 
being those due to the definition of the variables). 
Finally, then, to the order of approximation contem- 
plated in this paper, it will be correct to solve Eqs. (12)- 
(15) and to apply at » = 1 the boundary con- 
ditions (22). From the definition of f(R, w) and 
the relation (11), it follows that these are equivalent 


to 
fw" gD qwll | 
. : = 0) 
© a is ig SS 
uc? M.4|3;.. u?/U |? 
, = =—(y + 1) ——] om | 
& | I . B° i .- A 
5) The Irrotationality Conditions 


Applying the Tschaplygin transformation to the 
irrotationality conditions 


and evaluating them at the plane w = 0, 7, one obtains 
the following useful relations: 


Blu.’ ]o = — 2 - [we], (24) 
l—s 
I ~ 
Blu,]o = 2 [w, Jo (25) 
l—s 
s? 
Blu, ]o = 2A a 
(1 — s*) 
‘ » | 
— (D274 + 9 + w) +2; [w,," Jo 
(26) 
l+s 
[vy,™]> = —B - [u,™ Jo (27) 
l+s 
[7 Jo 2s [w,” | (28) 
sl 











332 JOURNAL OF THE 


‘ARROW-HEAD"’ SWEPT BEHIND MACH 


CONE 


EXAMPLE 


(1) Statement of Problem and Formulation of Boundary 


Conditions 


Fig. 2 shows the configuration of the body—an “‘ar- 
row-head”’ at zero angles of attack and yaw with respect 
and Mach Number 


to a free stream of velocity U 
M.. ¢€is the thickness parameter. 


After certain geometrical maneuvers, the boundary 


condition at the body for y > 0 may be written as: 


[ez]o = el¢r]n — € cot A[y,], (29) 


Assuming analyticity of solutions near the- plane of 


symmetry, one writes 


¥, 2) In = [e2]o + € ix, y) + 


eyo(x,y) +... (30) 


Substituting expressions (8) and (30) into Eq. (29) and 


[o.“”( (x. ¥ 


collecting terms of like order in e, 


le. Jo = (31) 


[¢: Jo = —-y + [¢z” Jo = COE Ale," lo (32) 

Eqs. (31) and (32) are applicable on the plane of sym- 
metry within the leading edges. 

Utilizing the Tschaplygin transformation, 

that s = 1 represents the Mach cone, and defining / as 


noting 


the s-coordinate of the leading edges, Eqs. (31) and (32 
may be written as 


fw] = U @w=0,s< Jd (33) 
[w?]o = —yils, w) + [wu ]o — 
1+7 
BS eM (= 0,5<) (34) 


By symmetry, 


[w™]> = (wO], = 0 w=0,1<s< 1) (35) 
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Fic. 2. ‘‘Arrow-head” wing swept behind tip Mach cone 
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Fic. 3. Diagram showing notation in Tschaplygin plane 


uw) is even in both s and y 
v™ is even in 2, odd in y 36 
w™ is odd in z, even in y 


On the circle s = 1, 


a 0% we), =O 


u* Mott uO /U > (97 
~~ ey sed lim a4 
i Pea see 


The differential equations applicable in the Tschaplygin 
plane are given by Eqs. (12)—(15). Use will be found 
for the irrotationality conditions (24) through (28). 

Fig. 5 represents the Tschaplygin plane, in which the 
problem will be solved. Stated in its simplest terms, 
the plan of attack will be as follows: 

(a) Find uw, v, and w", which satisfy differen 
tial Eq. (12), and the irrotationality conditions, subject 
to the boundary conditions developed above. 

(b) Find the particular integral for u‘’, using the 
expression (21). To this will be added the complemen 
tary function required to ensure satisfaction of the 
boundary conditions, thus completing the solution for 
a, 

(c) Using these results, apply Eqs. (4) and (5) to 


find the pressure coefficient. 


2) Singularities 


Certain difficulties are encountered in carrying out the 
above procedure because of the presence of singular 
' and v” and, hence, in all components in the 
second approximation, at the leading edges. This 1s 


due to the subsonic nature of the flow at these edges 


ties in u 


and the consequent stagnation in the component of 
velocity in the plane of symmetry and normal to the 
leading edge. A perturbation theory always fails t 
represent stagnation conditions and provides singular 
ties in velocities instead. 

One would expect that the singularities arising in this 
problem at the leading edges would be of the same types 
yawed 


as those appearing in the case of an infinite 


wedge, when the angle of yaw is sufficient to provide 
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subsonic normal flow at the leading edge. This nor- 
mal component would contain the only singularities. 
Further, the subsonic two-dimensional flow over a 
wedge will have the same type of singularity as that of 
the infinite yawed wedge; the incompressible (1/7. = 
() two-dimensional flow over a wedge will afford the 
same singularities as the compressible subsonic case, in- 
asmuch as the leading edge in the latter case represents 
a stagnation point in whose neighborhood the local 
Mach Number is nearly zero. 

Examination of the incompressible, two-dimensional 
flow over a wedge shows singularities at the leading 
edge of the logarithmic type in the first approximation 
and the square of the logarithm and the logarithm in the 
second approximation. These singularities are all 
integrable. These, then, are the singularities to be 
expected in the ‘“‘arrow-head" case. 

Admission of such singularities raises the question 
of uniqueness. It is known that a potential problem 
has a unique solution when suitable boundary condi- 
tions are prescribed on a closed contour, provided the 
solution is to be regular everywhere inside. If singu- 
larities are admitted, an ambiguity may be expected 
i.e., certain harmonic functions having singularities may 
be introduced without disturbing the boundary con- 
ditions. One therefore inquires what harmonic func- 
tions having singularities of the type In |s — / and/or 
(In |s — /|)* on the axis of symmetry may be added to 
the regular solution for uw or u™ satisfying the ap- 
propriate boundary conditions, without disturbing these 
boundary conditions. Thus, in accordance with Eqs. 
(24), (36), and (37), it is required that such a singular 
function have zero angular derivative on the axis of 
symmetry, be even in both y and z, and vanish on the 
unit circle. It may be shown that the only singular 
function meeting these requirements is the following: 


Ky [In (rz) + In (73) — In (dry) — In (Ir4)] (38) 


Similar reasoning leads to the following ambiguous 


term for?! orv” 
Ky [In (rz) — In (73) — In (ry) + In (Iry)] (39) 


The constants A,’ and A,” have superscripts indi- 
cating the order of approximation involved. 

Since w'' and w' have boundary conditions on 
value, they can have no ambiguity of this type. There- 
fore, knowing w"’ and w’, one may find K,", K,°’, 


K,, and Ko. 


this purpose. 


Eqs. (25) and (26) are convenient for 


It appears that the lack of uniqueness discussed 
above arises because of the abandonment of the ve- 
locity potential in order to take advantage of the sim- 
plicity afforded by the Tschaplygin transformation. 


3) First-Order Solution 


(a) Solution for w"’. The differential equation is 
given in Eq. (12), and the boundary conditions are given 
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in Eqs. (33), (35), and (37). By analogy with incom- 
pressible fluid theory, the solution may be written 
down immediately, using the stream function for the 
incompressible source : 

w)(0; + 0. + 0; + 6, — wr) (40) 


wo = —(U 


The differential equation is 
The boundary conditions are 


(6b) Solution for u". 
given in Eq. (12) above. 
those given in Eqs. (25), (36), and (37), and from Eq. 
(24): 


[u]lo = O (Ss < 1) (41) 


It is plain that the solution analytic inside the unit circle 
is identically zero, leaving only the constant K,‘” to be 
found. Applying Eq. (25) to Eqs. (38) and (40), and, 
for convenience, evaluating terms near the leading edge, 
one may immediately write 

. 
wu = - [In (v2) + In (73) — In (/r,) — In (lry)] 

Bril—F- 

(42) 


(c) Solution for v'’.-By a procedure similar to the 


one employed to find u”, it is found that 


gv) = (U/r)[((1 + 2)/(01 I?) | [In (r2) — 
In (73) — In (/r,;) + In (lry)] (48) 


(d) Information Obtained from the Above Solutions. 
Expanding (40), one finds that 


y, = —U(6/r)(A + P)/] (44) 


and from Eq. (42) is obtained 


; v/U 1 (1 ’ 1 
lim k t | aha ( . : aa pes w] (45) 
o>: —_ J pT “ 


Eq. (45) may be substituted into the results of Light- 

hill’ for the strength and location of the shock wave. 
For use in finding vu, one puts wu", v, and w" into 

expressions (20) 


complex form to conform to the 


above. 


(4) Second-Order Solution 


In philosophy, the solution for « is obtained in the 
same way as was “‘' in the preceding paragraph. In 
detail, the procedure is as follows: One may write 


u = uy + u, | 


(Ak 
wi?) = w, + w, 2 f 46) 
f 


is obtained from Eq. (21), w,‘ is obtained 
and w," satisfy 


where u,“° 
from a similar expression and u,“ 
Laplace's equation in the Tschaplygin plane. 
(a) Find u,“; note that the functions obtained 
from the integrations indicated in Eq. (21) may be 
complex. «,‘ may be found in real form by adding to 
these complex functions any convenient functions of ¢ 
or ¢ alone. 
(b) Find uw, such that wu," + u satisfies the 
boundary conditions (24), (36), and (37) and has no 
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singularities of order higher than (In |s — /|)*. Ex- 
pression (38) is to be included in 7,“’. 

(c) Find w,“. 

(d) Find w,° 
boundary conditions (84), (35), and (36) and has no 
s —l1|)*. Note 
that no boundary condition on w at the unit circle is 


such that w,” + w, satisfies the 


singularities of order higher than (In 


applied. 

(e) Apply Eq. (26) to find K,”. 
fied by considering only those terms that are singular 
The equation will contain terms 


This step is simpli- 


at the leading edge. 
singular at ¢ = /, and one may require the equation to 
be satisfied upon arbitrarily close approach to the point 
f=. 
be dropped, and all nonsingular factors in singular 
The term involving 


All terms nonsingular at this point may then 


terms may be evaluated at ¢ = /. 
the unknown K, has a derivative with a simple-pole 
type of singularity. Thus, one needs only to formulate 
Eq. (26) for terms having simple poles in derivative. 
Finding A,’ in this way makes it possible to ignore 
any boundary condition on w'” at the unit circle, since 
its application would result in a regular function that 
would make no contribution to the determination of 
K,. 

This procedure has been carried out and formulas 
have been obtained from which a special case has been 
computed. These formulas are extremely lengthy 
and are not particularly edifying in their complete form. 
They have therefore been omitted from this paper. 


Results appear in Fig. 4, where the free-stream Mach 
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Number is taken to be V 2; the free-stream velocity 
unity; and the ratio of specific heats, 1.4. The oie. 


foil geometry, in accordance with Fig. 2, is fixed py 


setting 2e equal to 6° and 2A equal to 57.6 
RESULTS AND DISCUSSION 


The results of a general nature applicable to a prob 
lem involving a body lying entirely within the tip Mach 
cone may be summarized as follows: 


(a) It has been shown how the results of Lighthil] 
may be adapted to the formulation of “boundary con 
ditions” on velocities at the Mach cone. 

(b) The differential equation for the second ap 
proximation has been formulated and its particular 
integral has been developed in integral form involving 
complex quantities in the Tschaplygin plane. 

(c) A method has been developed to express body 
boundary conditions at some mean plane, accurate to 
any order of approximation. 

(d) <A rather complete discussion is presented con 
cerning singularities arising due to the presence of a 
certain type of stagnation point. 


The techniques of this paper could be extended to 
treat the case of a body lying in part beyond the tip 
Mach cone (e.g., a rectangular wing tip). However, 
in such cases, a successive approximation approach of 
this type appears to lead to infinite shock strengths at 
the ‘triple points” arising in this type of flow.* Until 
this difficulty is surmounted, it would appear imprudent 
to attempt pressure calculations on bodies extending 
beyond the tip Mach cone. 

The results of the present method for the “arrow 
head” at zero yaw and lift--as they appear in Fig. 4 
show a rather small second-order correction to both per 
turbation velocity and pressure coefficient. This cor 
rection is in the opposite direction to that found for the 
two-dimensional wedge, due, perhaps, to the fact that 
the leading edge is supersonic in the latter case and sub 
sonic in the former. 

It may be noticed that the dominant second-order 
singularity asserts itself only on very close approach to 
the leading edge, a circumstance which augurs well for 
the reliability of the solution in the remainder of the 
flow field. 

The size of the second approximation computed for 
a rather thin “‘arrow-head”’ indicates that a certain im 
portance should be attached to second-order contribu- 
tions to drag, ete. For a limited class of problem, the 
results obtained from the present theory might syste- 
matically be computed to provide data of engineering 
usefulness. 

Although his expansion takes a form different from 
that used here, the rapid and accurate convergence of 
the cone solution obtained by Broderick* argues well 

* Lighthill,? pp. 1222 and 12238. 

(Continued on page 383) 
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The Incompressible Boundary Layer with 
Pressure Gradient and Suction 


LEON TRILLING* 
Northrop Aircraft, Inc. 


SUMMARY 


This paper presents a generalization of the Blasius series to a 
boundary layer with arbitrary pressure and suction distribution 
By transforming the boundary-layer equations from the physical 
x, yplane to the x, « plane, one formulates the problem as a single 
nonlinear parabolic partial differential equation with given initial 
and boundary conditions. This equation is solved approxi 
mately, and a method is shown for finding the stability 
threshold. Three the 


boundary layer, the Schlichting boundary layer, and the boundary 


new 


special examples are given: Blasius 


layer along the upper surface of a laminar flow airfoil at an angle 
of attack 
to give completely stable flow over the positive pressure gradicut 


In particular, a suction distribution is constructed 
region of the airfoil 


NOTATION 


Y, ) = space coordinates 

L length parameter (airfoil chord) 

l = free-stream velocity 

l = velocity component parallel to plate surface 

J velocity component normal to plate surface 

p density 

r dynamic pressure 

y kinematic viscosity coefficient 

R overall Reynolds Number, R = Ul L/, 

x dimensionless chordwise coordinate, x = x/L 

y dimensionless normal coordinate, y = (}/L) VR 

u = dimensionless velocity component parallel to plate 
surface, u = U/L; 

r dimensionless velocity component normal to plate 
surface, v = (_V/Uo) VR 

p dimensionless dynamic pressure, p = P/pl 

¢ = shear parallel to plate, ¢ = uy 

¢ = shear at the surface of the plate, go = uy(x, 0 

6 displacement thickness of the boundary layer 

7] = momentum thickness of the boundary layer 

D = friction drag of airfoil 

é = velocity at the inner friction layer 

Réer stability threshold 

R5 = local Reynolds Number based on displacement thick 
ness, Rs = ub V/R 

a st ibility index, a = Rs, f/Rs 


INTRODUCTION 


7 RAPID DEVELOPMENT of high-speed aircraft has 
revived the chronic interest of designers in de- 
creasing the friction drag of airfoils. One cause of the 
high friction drag is transition of the boundary layer to 
the turbulent régime well before the trailing edge. 

It is known! that suction of a small quantity of fluid 


Irom the boundary layer retards transition. There- 


Received August 31, 1949. 
* Consultant in Aerodynamics. Also, Instructor in Applied 
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fore, systematic investigations of boundary-layer flow 
with suction were undertaken independently in various 
countries. They include theoretical discussions of 
flow with various pressure and suction distributions* ° 
and experimental studies of boundary layers with suc- 
tion slots!’ !* and with continuous suction.'* Some of 
that research was described by Goldstein in the Eleventh 
Wright Brothers Lecture."' 

In much of that work, however, there is no attempt 
to determine the stability of the flow. 


families of velocity profiles are sometimes introduced, 


One-parameter 


and the boundary layer is analyzed by means of the 
von Karman-Pohlhausen method.'® Such a procedure 
does not lead to an accurate estimate of the stability of 
the flow, since it gives no precise information about the 
velocity profile on which the stability is critically de- 
pendent.!® 

A method to develop a generalized Blasius series is 
therefore presented below. In particular, the Blasius 
boundary layer! and the Schlichting asymptotic bound 
ary layer* are discussed, and a detailed study is made 
of the boundary layer along the upper surface of a 
laminar flow airfoil (NACA 66-2,015) at an angle of 
attack (C;, = 0.2) with various suction distributions 
in the region of positive pressure gradient. 

These calculations give the boundary-layer profiles at 
various points along the chord, from which a critical 
Reynolds Number is determined, below which the flow 
is completely stable. But experimental results'”'**6 
indicate that the flow remains laminar well beyond that 
Reynolds Number, probably because the disturbance 
components of unstable frequencies grow so slowly that 
the affected fluid particles travel some distance down- 
stream before the transition stage is reached. Thus, 
one cannot predict transition, but one can compare the 
stability of various flow patterns and rely on flight tests 
for the actual transition Reynolds Numbers. 
BOUNDARY-LAYER EQUATIONS AND THEIR 
APPROXIMATE SOLUTION 


(1) THE 


Consider the steady flow of a real incompressible 
fluid past a semi-infinite, infinitely thin, two-dimen- 
sional flat plate parallel to the direction of undisturbed 
flow. The pressure distribution is given along the sur- 
face of the plate, which is made of porous material so 
that some fluid can be sucked through it. 

In a system of Cartesian coordinates with the origin 
at the leading edge of the plate and the x axis along its 
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surface, the following equations describe the motion of 
the fluid in the thin layer adjacent to the plate: 


ur +v, = 0 (1.1la) 
UU, + Vy + pr = Uy (1.1b) 
py = 0 (1.te) 
with the boundary conditions 
u(x, 0) = 0, w(x, O) = a(x) (1.2a) 
v(x, ©) = u,(x, ~) = 0 (1.20) 
u?(x, ©) + 2p(x, ~) = 1 (1.2c) 
u(O, y) = 1 (1.2d) 


The coordinates and the physical parameters are made 
dimensionless in the conventional manner (see ‘‘Nota- 
tion’’). 

Since p is independent of y, Eqs. (1.1)—( 
tem of two equations to determine u and v. 
quires u, which determines the drag and stability prop- 
erties of the flow. 

A coordinate transformation due to Crocco 


.2) is a sys- 
One re- 


20 


is now 


introduced: 


(1.3a, b) 


~ 
= 
Il 
~ 
= 
—_ 
= 
< 
Str 
II 
= 


with the Jacobian 


J = O(u, §)/O(x, y) = u, (1.3c) 


The Jacobian is finite everywhere in the boundary 
layer except at the leading edge, where a discontinuity 
in u is needed to satisfy boundary conditions (1.2a) 
and (1.2d) simultaneously. It vanishes at the outer 
edge of the boundary layer and at points where the 
flow separates from the plate. Since in the boundary 
layer u(y) increases monotonically from 0 to u(x, ©), 
the transformation (1.3) is one-valued. The Jacobian 
vanishes as y becomes infinite because the boundary 
layer has no clearly defined outer edge. 
at any separation point where the flow breaks down. 


It also vanishes 
In terms of the new dependent variable g(u, £) = u,, 
the system (1.1)—(1.2) is transformed into the equation 
, ” 
uge — PP (E)eu = Pun (1.4) 
with the boundary conditions 


(1.5a) 
(1.5b) 


¢.(0, g) a [p’(&) ¢(0, £)] = v(E) 
y[u(é), £] = O 


To specify the problem completely, one must pre 
scribe an initial condition by defining the singularity 
at the leading edge; this was investigated by Carrier 
and Lin,”! who found 


lim (0, &) = 0.382/VYé (1.5¢) 


tE— 0 


Eq. (1.4) is a parabolic second-order nonlinear equa- 
tion that describes the diffusion of vorticity from the 
plate surface into the boundary layer. The charac- 
teristic lines are £ = constant, and the boundary condi- 
tions at £ 2 & have no effect on the solution in the re- 
gion — < &, since no signal can propagate upstream. 
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The flow along a finite plate is therefore identical] with 

the flow along a finite portion of a semi-infinite plate. 
No general method of exactly solving the systey 

(1.4)-(1.5) was found, but an approximate solution js 

developed here. The solution is expanded into a Mac 

laurin series: 


g(u, &) = DP On(§)u" (1.6) 
0 

which is assumed to converge in the interval 0 < 
u(é). If the series is then substituted into Eq. (1.4) 
and boundary condition (1.5a), one can determine all 
the coefficients ¢,(£) in terms of the shear at the plate 
surface go(£) and its derivatives by equating to zero 
the coefficients of all the powers of uv. The values of 
To determine 
¢o(£), one substitutes the series (1.6) into boundary con 


¢n(¢go) are tabulated in Appendix (A). 


dition (1.5b), with ¢, given in terms of g and with 
u(é), a known function of & defined by the potential 
flow. 

The order of the resulting equation for gp» increases 
with the number of terms taken, and, to solve the prob- 
lem rigorously, one must consider an equation of in 
finitely high order. The for that 
Actually, in the calculations 


initial condition 
equation is Eq. (1.5c). 
given below, six terms are taken, so that one has a 
second-order equation; condition (1.5c) is satisfied with 
a small numerical error. 

When ¢» and ¢, are known, the coordinate y is found 


“du 7 
Be = (1a) 
0 ¢g 


This integration is carried out by constructing the series 


from 


x 


l/e = > ®,(é)u" 


0 


(1.7a) 


convergent in the region 0 < u < u(£) but divergent 
at u = u(§) because of Eq. (1.5b). This series is inte 
grated term by term, and y is expressed as a power series 
in “, which diverges when u = u(&), since u takes the 
value u(é) when y is infinite (see reference 22). 

The velocity profile u(y) is computed by inverting the 
series y(u) (see reference 22). The velocity profile is 
then given by 

u(x, y) = pi [v,,(x) n!|y" (1.8) 
] 
The coefficients Y, are tabulated in Appendix (A). 

This method leads to results similar to those given by 
Green,” but it is expected that the application ol 
Crocco’s transformation allows a more accurate deter- 
mination of the shear ¢p. 

Other important properties of the boundary layer are 
its thickness, drag, and stability. These are now dis- 
cussed. 

The displacement and momentum thickness are de- 
fined by 
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6 = [1 /u(x)] f [u(x) — ul] dy 
0 


6= j1 mort f [u(x) — ujudy  (1.9b) 
0 


(1.9a) 


In the wu, plane, these integrals become 


| Rites du : 
§= [u(é) — ul (1.9a’) 
u(é) Jo 7) 
; | sf ‘te du 9b’ 
= A u(g) — ul u (1.9D ) 
[u(é) |? Jo ¢ 


The integrand is regular throughout the range of inte- 
gration, since (see references 23 and 24) 


im u/uy = —2xu(x)/y = 0 (1.9¢) 

The boundary-layer thickness is, therefore, 
6 = Ey —slu(x)) n(n + 1) (1.10a) 
= ob, _s[u(é)]" (n + 1)(n + 2) (1.10b) 


The coefficients ®, are tabulated in Appendix (A). 
The drag is given by the integral 


" 
D(x) - f go  E)dé 


The stability of the boundary layer is discussed in 
terms of a stability threshold R,., and a stability index 
a. The stability threshold is defined as the largest 


(1.11) 


Substituted into Eq. (1.12), these give 


~ ®,(&) m 2 
rel De +3 é oot! — ac|| S (2 + 1) Ons er] =: (0.58 b ele" 
0 


9 n+ 1 


0 
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Reynolds Number based on displacement thickness at 
which any small disturbance of the laminar flow is 
damped out. Lin'® has related the instability of a vis- 
cous flow to the existence of an inner friction layer and 
found the dependence of the stability threshold on the 
velocity c of that layer. He has also determined the 
velocity c(x) in terms of the velocity profile u(x, y) of a 
boundary layer. The series (1.6) is a particularly con- 
venient tool to investigate the stability threshold of vari- 
ous flows. 

The velocity c of the inner friction layer satisfies the 
following equation :*” 


to[(2goyv Gg=- 3] (cCyy Cy”) = 0.58 (1.12) 


and the stability threshold is approximately equal to 


R;., = 25gou*(x) d(x) /c4(x) (1.13) 


Eq. (1.13) differs from Eq. (2.27) of reference 27 by the 
factor u‘6 because of differences in the definitions of the 
normal dimensionless coordinate and the velocity u. 

The following expansions are useful in solving Eq. 
(1.12) approximately : 


“du — #,(é) 
vio) = f — = >: é er (1. 14a) 
‘ 0 ¢ 90 at i 


x 


or, replacing ¢, and ®, by their values and ordering in increasing powers of c, 


0.1845 go? + 1.5535 go(voeo + p’) ¢ + 1.2768 


Yo 


(0.6578 Voeo Yo) + (0.0602 p’ go’ ¢go~) 


Voyo f ’ aT j . , ” 
— 2 P’ (1.1845 ve2ge? + 0.1667 oop’ — 0.4255 p | c? + 11.0461 [v0’ + (p"/0)] — 


ty - g(c, g) = ; I ¢n(§)c" (1.14b) 
0 
Cyy = o¢u(c, §) = b ee || S (n+1)¢n4 ee] 
0 0 
(1.14e) 
(1.15a) 
(Vogo + p’)(2v0g0 + p’) Cc + | 0.5922 Yo’ + 
$0 
4 Voyo + p (—0.1667 Uo? oo* + 0.2872 Vo7go7p’ + 
Yo 
1.294 vogop’? + .9025 p') hes +.-.-=0 (1.15b) 


Since, in most cases, c < 0.4, a good approximation is obtained by keeping only four terms in the series (1.15). 


In this manner, an approximate value of c(x) is ob- 
tained by solving a fourth-degree algebraic equation, 
whose coefficients depend on the parameters of the prob- 
lem p’(), vo(€) and the shear at the wall g. This 
value of c is then substituted into Eq. (1.13), and the 
stability threshold R,., is computed. 

To determine whether a given flow is completely 
stable, one must compare its local Reynolds Number to 


its stability threshold; the local Reynolds Number of a 
boundary layer based on the displacement thickness is 


R; = u(t)5(E) VR (1.16) 


where R is the overall Reynolds Number of the flow 
R = UL/». 


If the stability index a = R,,,/R, is greater than 1, 


the flow is completely stable; if it is less than I, the 
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flow is unstable for a certain range of disturbance fre- 
quencies. The rate of growth of probable frequencies 
must then be investigated before the transition point 
can be predicted. Little is known about that problem, 
but a comparison of stability indexes gives a qualitative 
idea of the stability of various flows. 


(2) THE BLASIUS AND SCHLICHTING BOUNDARY LAYERS 


To illustrate the method outlined above, consider the 
boundary layer along a flat plate without pressure gra- 


dient or suction. The equation of motion is 


Ud: = O° Pun (2.1) 


with the boundary conditions 


g,(0, £) = 0, g(1) = 0 (2.2a, b) 


and the initial condition 


lim ¢(0, £) = 0.332/+/¢ (2.2c) 


t— 0 
Crocco” shows how this equation can be solved by 
separation of variables to give the familiar Blasius 
series. Below, the series (1.6) is constructed, and the 
following equation is found for ¢: 


go’ go" £0 - 4g"? 
£0 a : 
6g" LSO Yo” 
go’ go" — 24 go" ¢0'¢o0 + 230" _ 23) 
12,960¢0° -_ 
The solution of that equation is 
go = 0.320/VE (2.4) 
The velocity profile is given by the series 
9.399 2. — (0-320) * 
u = 0.62 = ° 
Vx 4! <* 
1] _¥ P 
—(0.320)? = +... (2.5) 
(! s* 


This is the Blasius solution with an error of 3.5 per cent 
The displacement thick- 
The sta- 


in the numerical constant. 
ness of this boundary layer is 6 = 1.73//E. 
bility parameters c, R;,,, and @ are 
R;., = 511, a = 281/V RE 


(2.6a, b, c) 


c = 0.409, 


The velocity c and the stability threshold R;,, are inde- 
pendent of the chordwise coordinate £ because of the 
similarity properties of the Blasius profile; the index 
a decreases as 1/+/£ because the boundary layer thick- 
ness increases. The value of the stability threshold 
agrees closely with that obtained by Lin (502) and by 
Chiarulli and Freeman (465). 

Schlichting* investigated the asymptotic behavior 
of a boundary layer with no pressure gradient and uni- 
form suction v, at large distances from the leading edge. 
All parameters are then independent of £, and the solu- 
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tion of Eq. (1.4) with boundary conditions (1.5) js 
Yo = (4 — 1 )v% (2.7 
so that 
y = (1/m) log (1 — 4), u=1l—e” (9¢ 


rhe same result is obtained by constructing the series 
(1.6). Since all the terms except ¢g» and ¢; vanish, the 
solution becomes 
u = >) v(y"/n!) = 1 — e” (29 
I 


= —1/v, and the sta- 


~ 


The displacement thickness is 
bility parameters are 
R;.- = 41,000, a@ = 41,000V,/U, 


(2.10 


c = 0.1565, 


This value of R;,, agrees well with the value 41,770 
found by Chiarulli and Freeman.*” Note that the 
velocities of the inner friction layer and the stability 
threshold are independent of the suction velocity. 
Since the boundary-layer thickness is inversely propor- 
tional to the suction velocity, the stability index is pro- 
portional to it. Thus, for any given free-stream ve 
locity, there is a critical suction velocity above which 
the boundary layer is completely stable. If the free 
stream velocity is 600 ft. per sec., the suction velocity 
is 0.0146 ft. per sec. 


(3) BouNDARY ALONG UPPER SURFACE OF AN AIRFOIL 


The method developed above is now applied to the 
investigation of the boundary layer along the upper 
surface of an NACA 66-2,015 laminar airfoil at a lift 
coefficient C, = 0.2 and a Reynolds Number R = 
36,000,000. 

The flow breaks down shortly after the pressure 
minimum if there is no suction. It remains laminar 
up to 80 per cent of the chord if just enough suction is 
applied to make u,,(0) slightly negative at the pressure 
minimum; the drag increases considerably if the 
boundary layer is overstabilized by too much suction. 
A suction distribution that maintains a stable low drag 
boundary layer over the rear portion of the airfoil is also 
constructed. In all this work, the airfoil is approxi- 
mated by a flat plate. 

The pressure distribution given in reference 28 is 


shown on Fig. 1. It is closely approximated by the 


expressions 
p’(é) = —0.022 (0 < — < 0.63) (3.1a 
p(t) = 0.751 (0.63 < — < 1) (3.1b 


Since the pressure gradient over the front portion ol 
the airfoil is small, the boundary-layer flow there is ob 
tained as a perturbation of the Blasius flow. The val- 
ues found for go (0.63) and go’ (0.63) are then used as 
initial conditions in computing the flow over the rear 


of the airfoil. By neglecting the discontinuities intro- 
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duced into the parameters p’ and v», one introduces dis- 
continuities into 6 and 6, but these turn out to be 1m- 
portant only in the immediate vicinity of the pressure 
minimum. 

The velocity of the free stream along the front portion 
of the airfoil is not unity but 1.20 because of the pres- 
sure distribution. It follows from dimensional consider- 
ations that the appropriate perturbation series in that 


case takes the form 


p’ li 
go = Uy’ | op + TNE) fo(d) + 
Uo” Uo! 
(3.2) 
where ov, represents the shear of Blasius flow. This 


series is substituted into the expressions for ¢,,(¢o), and 
into the condition (1.5b). If six terms are taken and 


the coefficients of the various powers of p’ are made to 
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vanish, it is found that ¢, satisfies the approximate 
Blasius equation, and the perturbation functions 
fn(§) all satisfy the following linear second-order equa- 
tion: 


0.0055&*f,," + 0.0385éf,’ + 
0.0375 fn + Ant” 


(3.3a) 


with the initial conditions 
lim f,(é) = o(¢ 2) dim fn (&) = of& va 


E— 0 E— 0 


(3.3b) 


The constants A, are evaluated, and the solution of the 
perturbation problem is found to be 


go = (0.414/1/E) + 0.025 VE + 0.00015 2°" 
(3.4) 


The shear, the displacement thickness, the stability 
threshold, the Reynolds Number based on displacement 
thickness, and the stability index are shown on Fig. 2. 
At the point of pressure discontinuity, one finds 


go = 0.525, go’ = —0.387 (3.5a, b) 


While the boundary layer is completely stable for a 
short * distance, flight tests show that, under normal 
conditions, a caretully finished surface supports laminar 
flow up to the pressure discontinuity’ ' or a dis- 
placement thickness Reynolds Number of order 6,500. 
Therefore, no suction is applied to this portion of the 


surface. 
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Fic.3 Flow ovER REAR OF AIRFOIL 


D with Constant SucTION Uo=-!.53 
2) with Uyy = - 0.15p' 








The flow over the rear portion of the airfoil is in- 
vestigated under several conditions. First, a constant 
suction velocity just sufficient to make u,,(0) slightly 
negative at £ = 0.63 is applied. Next, the amount of 
suction is doubled. Finally, a suction distribution 
maintaining a laminar boundary layer of moderate 
drag is constructed. 

The pressure gradient is too large for any perturba- 
tion procedure to be practical. Instead, the solution is 
developed into a Taylor series about the point — = 0.63. 
When the suction velocity v is — 1.53, the Taylor series 
is 
go = 0.525 — 0.387 (€ — 0.63) + 


0.166 (& — 0.63)? ++... (3.6) 


A second expansion of the solution about £ = 0.83 is of 
the form 


yo = 0.455 — 0.321(— — 0.83) — 


7.85(— — 0.83)? +... (3.6a) 


indicating that the flow is breaking down, since go — 0. 
This boundary layer is illustrated on Fig. 3. 

Since the velocity v) = — 1.53 is not sufficient to main- 
tain a laminar boundary layer over the rear portion of 
the airfoil, an appreciably higher suction velocity is in- 


The resulting equation is 


{ , — Apis , Apis? 


: go" 2¢0' 608 


A[4A* — (29/3)A? + (11/3)A — (1/6) pu’ — A[30 A* — 111A4* + 824? — 134 + (1/4) ]p/%u! 


20g!” 


2 
¥0 


The substitution 


go? = up’f(log u), 


ou? (1 3A — (1/2) |p’u 

¢o't {" 4 [ iv lp 4 
6 12g? 

[42 A* — (287/2)A? + (71/4)A — (17/4)] pu) 

180 go® f 


vestigated. When 7 = —3.00, the Taylor series be- 
comes 
go = 0.585 — 0.387(— — 0.63) + 245(E — 0.63): 


(3.7 
The shear increases rapidly, and the boundary layer 
becomes extremely thin. While the flow remains lami- 
nar to the trailing edge, the profile drag becomes ey. 
tremely large. 

It is therefore clear that the suction must be adjusted 
to local conditions so that the boundary layer will re. 
main laminar without undue drag increase. One way 
to obtain such a result consists in keeping the curvature 
at the wall ,,(0) slightly negative and constant. The 
suction velocity and the shear are then connected by the 
relation 


(3.8 


goo + p’ iad —Ap’ 


where A is a small positive constant. A differential 
equation for go is constructed by substituting Ey, 


(3.8) into the parameters ¢, of Appendix A and by using 


for the velocity at the outer edge of the boundary layer 
the Bernoulli relation, 


/ 9 ‘ ; ‘ 
u = Vu? — 2p't (3.9) 


A(2A — 1)p’*us  A[5A* — 44 + (1/2) ]p/*ut . 


12¢g0° 
st { 
180¢0!? ct 
[5A? — (29/6)A + (1 6) |p’*u? + 
20 go* 
gogo” aa 4(¢o')? 


> = = 0 
¥ 180 go 


(3.10 


g = df/d(log u) (3.11) 


transforms Eq. (3.10) into the first-order total differential equation 


124 — 1 
12f 


45A + 2 50A?2 — 114A — 1 
-120f? — 240f 3 
A(24A* — 58A? + 22A — 1) _ 
120° . 
30A?2 — 29A + 1 
oz SC 


i 

' Ar 4 90A ae 168A* — 

of \6 360f 
with the boundary condition 

g(0.310) = 0.525 (3.13) 


followed by the quadrature 


Reins -f. df 
” 1.2 Josio g 


Some information on the nature of the flow is obtained 


(3.14) 


_ A(10A? — 8A + 1) 


24ft 
A(120A* — 444A? + 328 A* — 52A + 1)) 
720f° f 
574A? + “1A — 17) _ (i >) =-0 (312 
720f* f—_-180f* 


df of 


by considering the solution when f — 0. The asymp- 
totic formula valid in that case is 


lim g(f) = 


f—mo 
[4 120A* — 444A* + 328A" — 524 +) 
6f? 
(3.15) 
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INCOMPRESSIBLE 


When A < 0, the quantity under the radical is positive, 
so that g is real; the flow may therefore break down, as 
it does when A is —1. When A > 0, g isimaginary, 
so that the flow does not break down. 
Numerical integration of Eqs. (3.12)—(3.14) was car- 
ried out for A = 0.15. The results are shown on Fig. 3. 
The velocity at the inner friction layer is given by 


0.1845 — 1.553549 + 1.2768A(A + 1)? — 
(0.2961f (f + g) + A(1.1845A? + 2.2023A — 
0.5922) In? — [0.8821f(f + g) + A(0.1667A* + 


0.787242 — 0.21944 + 0.0624) |n? +... = 0 (3.16) 


where » = p’c/go”. Note that, when A = 0.15, there 
is no inner friction layer, so that the flow is stable at all 
Reynolds Numbers. An examination of Eq. (3.16) indi- 
cates that this result is due to the fact that the function 
gis not positive, so that the shear does not decrease as 
the fluid flows downstream. 


CONCLUSION 


In the preceding pages, an attempt was made to 
work out a method of determining the stability of a 
boundary layer with arbitrary pressure gradient and 
suction distribution. By transforming the equations 
of motion from the x, y plane into the x, u plane, one 
can reduce the problem of the boundary layer to the 
solution of a single nonlinear parabolic equation for the 
diffusion of the shear. That equation is solved ap- 
proximately by means of a Maclaurin series expansion, 
and the resulting solution is especially convenient for 
stability calculations. 

In particular, it is shown that along the rear portion 
of a laminar flow airfoil, in the presence of an adverse 
pressure gradient, it is possible to keep the boundary 
layer completely stable by the proper suction dis- 
tribution. This turns out to be such that the curvature 
of the velocity profile at the wall is kept slightly nega- 
tive. While the drag is greater by 8 per cent than the 
equivalent Blasius drag, it does not decrease toward 
the trailing edge, so that this distribution may not be 
suitable for large proportions of the chord. 
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A Study of the Nonlinear Characteristics of 
Compressible Flow Equations by Means 
of Var ‘ational Methods 


CHI-TEH WANG* 


anpG.V. R. RAOF 


New York Unwersity 


SUMMARY 


By means of variational methods, the problem of a steady 


jrrotational compressible flow past an arbitrary body may be 
solved directly in the physical plane. The resulting equations 
we nonlinear in nature. A study of these equations is made in 
this paper. It is found that at high Mach Numbers the solu- 
tion may no longer be unique. Beyond a certain limiting Mach 
Number, no physically possible solution exists. At this limiting 
yalue, it is shown that an unsymmetrical flow pattern may ap- 


pear 
It is also shown that, if only the velocity and pressure distribu 


tions on a given profile at various Mach Numbers are to be 
found, good approximate results may be obtained by the use of 
the linearized variational integral. Comparisons are made be- 
tween the linearized and nonlinear solutions in the case of a circu 
lar cylinder, a Kaplan’s bump, an elliptical cylinder, and a 
It is seen that in all these cases the linearized solutions 


sphere 
By working with the 


give close agreement to the nonlinear ones. 
linearized integral, the labor saved is enormous. 


INTRODUCTION 


[' PREVIOUS PAPERS,'~* one of the authors has applied 
the variational method to compressible fluid flow 
problems following the Rayleigh-Ritz procedure. The 
results show surprisingly good agreement between the 
variational method and the other approximate meth- 
ods. A study of the results also reveals that the vari- 
ational method gives good approximate solutions at 
both high and low Mach Numbers and for flows past 
both thick and thin bodies. 

In carrying out the calculation, 
linearization of the fundamental equations is neces- 
sary; also, the calculation is made directly in the physi- 
cal plane. It therefore gives the stipulation that a 
study of such solutions may throw some light on the 
nonlinear characteristic of the fundamental equations. 
An investigation in this direction has indeed proved to 
be fruitful. At high Mach Numbers, it was found that 
the solution of the equations may cease to be unique. 
Mach Number, no physically possible 
solution seems to exist. At this limiting Mach Num- 
ber, the occurrence of an unsymmetrical flow pattern 


it is seen that no 


Beyond some 


is shown to be possible. 
The study reveals that, up to this limiting 
Mach Number, if the variational integral is linearized 


also 
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by assuming 
¢= o> + ¢ and ¢& < ¢,; 


where ¢ is the velocity potential of the flow, ¢, is that 
for the corresponding incompressible flow, and ¢, is the 
correction part due to the compressibility effect, the 
linearized solution gives good approximation to the 
nonlinear one. This is shown by comparing the results 
obtained by the linearized and the nonlinear integrals in 
the case of flow past a circular cylinder, a Kaplan’s 
bump, an elliptical cylinder, and a sphere. 

METHOD 


THE VARIATIONAL 


By denoting the velocity vector by q; the Cartesian 
coordinates by x), %2, and x3; and the corresponding 
velocity components by 7, v2, and v3; the motion of a 
steady irrotational flow is to be determined from the 


following fundamental equations: 


Continuity equation div (eg) = 0 (1) 
Condition of 
irrotationality rot ¢ = 0 (2) 
Equation of motion p 0q°/Ox; = —20p/Ox; 
(i = 1,2,3) (3) 
Equation of state p = A + Bp’ (4) 


v2” + v3"; pis the pressure; p is the 


and a repetition of 


where gq? = 2" + 
density; A, B, and y are constants; 
subscripts indicates summation. 
It can be shown? that the first variation of the in- 
tegral 
S pdr (5) 
gives Eqs. (1), (2), and (3) as its Euler’s equations 
in the Calculus of Variation if Eq. (4) is given and if 


¢ poog,dS = 0 (6) 


In integral (5), dr is an elementary volume, and the 
integration is to be extended over the whole volume of 
the fluid. In integral (6), dS is an elementary surface 
on the boundary of the fluid, g, is the component of the 
velocity normal to the boundary, 4¢ is the variation of 
the velocity potential ¢, and the integration is taken 
over the total boundary surface. If integral (6) is not 
zero, suitable integrals must be subtracted from Eq. 
(5).>* 
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In the two-dimensional problem of a compressible 
flow past circular cylinders, ¢ may be assumed to be in 
the following form: 


1 , 
in u(r + ) cose + 6+ 
r o 


Tv 
l l 
mr” (m+ 2)r™t2 x 
m n 


(A nm, cos nO + By, sin n@) (7) 


where U is the velocity of undisturbed flow; «x is the 
strength of circulation; A,,,, and B,,, are undetermined 
parameters; and 7 and @ are the polar coordinates. It 
may be noted that Eq. (7) satisfies the boundary con- 
ditions—i.e., atr = ~,g = Uandatr = 1,0¢/or = 0. 

For any arbitrary shape, one may map the circular 
cylinder in the complex z-plane into the desired shape 
in ¢-plane by the conformal mapping function ¢ = 
f(z). It may be proved! that $(¢) satisfies the boundary 
conditions. 

With ¢ as assumed in Eq. (7), the variational integral 
(=* nav he put into the following form: 
I S° lam? — 9(2)?2/|de/d2|2]”~ ? | dg/dz|? X 

s-plane 


r dr d0 + 2yx(Gm? — U?2)"“O-YDUAy/(y — 1) (8) 


where the last term is due to the fact that integral (6) 
is not zero and 


da’ — U? + [2a0?, (y = 1)| (9) 
q(z)? = (O¢/Or)? + (O¢/rd0)? (10) 
g(f)? = g(z)?|dz/d¢|? (11) 


dy being the velocity of sound at infinity and g,, being 
the maximum possible velocity of the flow which fol- 
lows from the fact that p cannot be negative. 

When a compressible flow problem is studied by solv- 
ing the differential equations, the condition of non- 
negative p is assured when the integrated form of 
equations of motion [Eq. (3)] is used. In the solution 
by the variational method, this condition requires that 
the integrand of / must be a positive definite function, 
which is equivalent to saying that g? must be less than 
dn” everywhere in the flow. Since this condition is not 
used in carrying out the computation, it must be in- 
corporated into the problem as an auxiliary condition. 


SOLUTION OF THE NONLINEAR EQUATIONS 


By following the Rayleigh-Ritz procedure, the vari- 
ational problem is transformed into an ordinary sta- 
tionary value problem. The solution of the problem 
then reduces to the solution of the resulting nonlinear 
simultaneous algebraic equations. 

The general solution of a system of nonlinear equa- 
tions with many real roots is a difficult problem. There- 
fore, only simple cases will be studied in this paper. 
Let us consider the flow past a circular cylinder with- 


out circulation. Then the flow must be symmetrica| 
Taking only one parameter, Ay, and using y = 2, th 
equation 0/J/0A;, = 0 gives the following cubic equa- 
tion: 


1.05556U* — (0.51852q,? — 1.82963U2)Ay + 
0.41111UAu? + 0.0507944,;* = 0 (19 


This equation can be solved by standard algebrai 
methods. Solving this equation, it is found that A, 
has three real roots when Mo < 0.615 and has only on 
real root when Mp) > 0.615. 

When My > 0.615, the one real root of Aj gives 
maximum velocities well above the maximum possible 
velocity gq», at all Mach Numbers and, therefore, is to be 
discarded. It indicates that when My > 0.615, there 
does not exist any physically possible irrotational 
flow. When My < 0.615, one of the three roots again 
gives maximum velocities well above q,, at all Mach 
Numbers and is again to be discarded. The magnitudes 
of the maximum velocities given by the other two roots 
are plotted in Fig. 1 against the Mach Number of the 
undisturbed flow. It is shown in reference 2 that the 
use of y = 2 instead of 1.405 has only slight effect on the 
magnitudes of velocities, but it is obvious from Eq. (9 
that the value of y has an appreciable effect ong, Ii 
the values of g,, corresponding to y = 2 are used, then 
only one of the two roots gives maximum velocities 
below q,. But if the values of g,, corresponding to 
y = 1.405 are used, it is seen that there exist two pos- 
sible velocities in the region from My = 0.537 to My = 
0.615. This indicates that in this region the solution 
may not be unique. 

Now let us try to solve Eq. (12) by the method oi 
successive approximation. Rewrite Eq. (12) as fol 
lows: 


Ay = (1.05556U? + 0.41111 UAy? + 0.050794A1)') + 
(1.03704a9? — 1.31111U%) (18 


where the relation g,,2 = 2a) + U* has been used. 

The method of successive approximation may be 
carried out as follows: First, assuming the nonlinear 
terms in A), to be zero, we have 


Ay, = 1.05556.M*/ (1.03704 1.31111 Mo") 


Next, substitute this value of A, into the right-hand 
side of Eq. (13). The value of Ay so obtained is then 
the second approximation. Repeat the process until 
the same value of A; is obtained in two consecutive 
cycles. This method works well and the rapidity o! 
convergence is satisfactory. Beyond M, = 0.61), 
however, the method fails. The reason for its failure 
now becomes obvious, because the only root of Aun i84 
large value and the first approximation is not a good 
approximation at all. 


It is usually believed that, once the symmetrical | 


flow pattern fails, the unsymmetrical flow patter! 
appears. By including an unsymmetrical term ™ 
¢, we have 
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The conditions 0//0A,;, = 0 and 0//OAy = O then give, 
respectively, the following two equations: 


1.05556 U* — (0.51852qm? — 1.82963U7)Aun + 
0.41111UAy? + 1.09630UA ys? + 0.050794A 11? + 
0.367464 yA 12” = @ (15) 

Ay|—(0.87037¢n2 — 2.67963U2) + 1.096830UAn + 
0.19021A 4,” + 0.35225A 12?| = () (16) 


An inspection of Eq. (16) shows that Ay, = Ois a root. 
With Ay, = 0, Eq. (15) reduces to Eq. (12), the solu- 
tions of which have been obtained previously. Elimi- 
nation of A,» between Eqs. (15) and (16) results in an- 
other cubic equation that can again be solved by the 
conventional method. The maximum 
responding to these roots, when Ay» + 0, are shown 
by the dotted curve in Fig. 1. It is seen that the curve 
Ay ¥ 0 meets the curve Ay» = 0 just about at the Mach 
Number where the latter curve turns. Physically, 
it indicates that, at this Mach Number, the unsym- 
metrical flow pattern is a possible solution. The actual 
occurrence of this flow pattern, however, must be de- 


velocities cor- 


termined from a consideration of the stability of the 
flow. 

From the above-mentioned facts, it is believed that 
the first occurrence of two solutions probably does not 
have much physical significance, because as Mo gradu- 
ally increases, the maximum velocity in the flow will, in 
general, increase continuously instead of at an abrupt 
jump. Also, the solution by numerical methods of a 
nonlinear differential equation may yield solutions 
other than the true ones of the original differential 
equation. The turning of the velocity curve and the 
appearance of unsymmetrical flow pattern do seem to 
indicate the failure of potential flow and the possible 
occurrence of shock waves, since, beyond this point, no 
physically possible flow exists. The value of My = 
0.615 is obviously much too high. 
the calculation only one parameter is used, and, there- 
fore, the results are accurate only to the degree that 
When six parameters are taken, it was 


This is because in 


order permits. 
reported in reference 2 that the failure of the method of 
successive approximation can be estimated to occur at 
about Mo = 0.46. 
the value estimated by other authors from other ap- 


This value is in close agreement with 


proximate methods. 

In the previous discussion, the values of qg,, corre- 
1.405 must be used if the discussion 
In carrying out the vari- 


sponding to y = 
is to have any significance. 
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ational method according to the Rayleigh-Ritz pro- 
cedure, y/(y — 1) must be an integer, and y was there- 
fore taken as 2. However, if the variational method 
is carried out according to the Galerkin’s procedure, 
there is no restriction on the value of y, and the correct 
value of 1.405 may be taken. Details of the applica- 
tion of Galerkin’s method to compressible fluid flow 
problems may be found in the papers by Wang and 
Brodsky.*'! By Galerkin’s procedure, the conditions 
o1/0Ai, = 0 and O//0A\. = 0 lead, respectively, to the 
following two equations: 


0.52778U — (0.51852a9? — 0.68200U?)Ay, + 
0.15781 UA)? + 0.58726L 'A 12” + 0.01 1440A,,3 + 


0.085770A A 12? = 0 (17) 


Ajo[— (1.57407ao? — 1.59391 U?) + 0.84544UAy, + 


0.12600A 1:7 + 0.24380A 127] = 0 (18) 


Solution of Eqs. (17) and (18) can be carried out 
similarly to that of Eqs. (15) and (16). The results are 
also plotted in Fig. 1. The lower branch of the curve 
obtained by the Galerkin’s procedure falls directly on 
the curve obtained by the Rayleigh-Ritz procedure. 
The other parts of the curves are not in such good agree- 
ment, but they all have the same characteristics and, 
therefore, the same discussion can be applied. This 
confirms the fact that y has only a slight effect on the 
resulting velocities, and the value of g,, corresponding 
to the correct y may be used in the discussion. It 
may be mentioned that the application of Galerkin’s 
method to such study has also been made by Akita,° 
but it appears that the method was not correctly 
carried out because of a mistake in the integration 


process. 


SOLUTION OF THE LINEARIZED EQUATIONS 


The assumed form of ¢ in Eq. (7) may be interpreted 
to consist of two parts—namely, ¢,, which is the ve- 
locity potential of the corresponding incompressible 
flow, and @,., which gives the correction due to com- 
pressibility effect. In the variational integral /, if 
we neglect terms of A,,, and B,,,, with higher than the 
second power, the final equations 


ol ol 
= ), = 
OA mn Oo B mn 


0, (et 0 = 1,2...) 


are then linear in A,,, and B,,,... We shall call the solu- 
tions so obtained the linearized solution. For example, 


the linearized solution of Eq. (12) is then 


Ay, = 1.05556.Mo"/ (1.03704 — 1.31111M]o?) 


An examination of the integral J shows that this 
process is the same as neglecting terms of ¢,, which 
has a power higher than 2 in the variational integral. 
The Euler’s equation of the linearized integral can be 
shown to be the same as the linearized form of the fol- 
lowing fundamental compressible flow equation: 


AERONAUTICAL 


SCIENCES—JUNE, 1950 
TABLE | 
Velocity Distribution of a Circular Cylinder at Jf, = 04 
, Rayleigh- 
g/4 Linearized Nonlinear Janzen 
/ U Solution Solution Method 
0 0 O 0 
10 0.3202 0.3194 0.3104 
20 0.6366 0.6349 0.6439 
30 0.9486 0.9410 0.9588 
40 1.2572 1.2470 1.2658 
50 1.5599 1.5522 1.5613 
60 1.8436 1.8450 1.8355 
70 2.0826 2.0974 2.0752 
80 2.2447 2.2712 2.2269 
90 2.3023 2.3335 2.2840 
TABLE 2 
Velocity Distribution g/U on the Surface of a Kaplan Bump at 
Mo = 0.83 
Perturba- 
x/ q Linearized Nonlinear tion 
U Solution Solution Method 
0 1.169 1.187 1.174 
0.1 1.163 1.180 1.168 
0.2 1.147 1.161 1.151 
0.3 1.116 1.131 1.245 
0.4 1.089 1.093 1.091 
0.5 1.052 1.051 1.054 
0.6 1.1014 1.009 1.017 
0.7 0.9784 0.9705 0.981 
0.8 0.9469 0.9401 0.949 
0.9 0.9226 0.9209 0.917 
0.975 0.9115 0.9180 0.892 
1.00 re a 0.882 


L(y — 1)/2] [4m? — (p/Ox;) (0p/Ox;)] — 
(Op/Ox;)?} 07p/Ox,Ox; — (Op/Ox;) (OG/OxX,) X 
(07¢/0x,0x,) = 0 (1,7, Rk = 1, 2,3) (19 


by assuming @ = ¢; + ¢, and ¢, to be small in compari 
son with ¢,;. The corresponding variational solution 
is then the approximate solution of the linearized equa 
tion. 
Let us again consider the case of flow past a circular 
Taking six parameters 
we have 


l l l 
u(- + ) cos + Ay ( ~ ) cos 6 + 
r r 3r 
l l } l l 
Ais _ 373 cos 36 + Az 3,3 — =,5 cos 6 + 
l l l 7 
-— — ]} cos 30 + Ais — — ]} cos 56+ 
. or” r 3r? 
( 3 | ) 
Asi -——J]cosé (20 
or dr’ 


0.4, the velocity distributions given by the 


cylinder without circulation. 


at M = 


linearized and nonlinear solution, together with thal 
obtained by Imai® using Rayleigh-Janzen method, at 


listed in Table 1. The coordinates of @ and r ar 
shown in Fig. 2. 
The circular cylinder may be mapped on to a sym 


metrical profile by the following mapping function: 


¢=2+ [(1 — d*)/2] + (d?/32?) (2 
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COMPRESSIBLE 


where d is a constant. This symmetrical profile is 
originally given by Kaplan’ and is referred to as Kap- 
jan's bump. Again, by assuming ¢ as in Eq. (20) and 
taking d? = 0.075, the velocity distributions obtained 
by the various methods are compared in Table 2. The 
value d? = 0.075 gives a bump with a maximum thick- 
ness of 5.13 per cent of the chord. The shape of the 
bump is shown in Fig. 2. The detail of the calculation 
of the variational method is carried out in reference 2. 

Asa third example, the flow past an elliptical cylinder*® 
is considered. Denoting the elliptical coordinates by 
tand 7, the velocity potential ¢ for the flow past the 
elliptical cylinder £ = & may be assumed to be 


> = gi + de 


o; = U(a + b) cosh (£ — &) cos n 
6. = [Au cos n/cosh (€ — &)| + [A1s3 cos* n 
cosh (¢ — &)] + [Asi eos n/cosh* (& — £&)| + 
[A33 cos* n cosh? (é — fo) | (22) 
where a and 6 are the semiaxis of the ellipse. With 


Eq. (22), the variational integral can be shown to be 


l= S'S (qm? — [(0G/08)* + (06/0n)*| 


1) 


c? (cosh? £ — cos? n)}" c? (cosh? & — 
cos? n)dé dyn + 2yrU ab(q,? — U2)“7-” x 
[An + (3/4)Ais3] /(y — 1) 


(23) 


where c? = a® — 6°. The velocity distributions on an 
elliptical cylinder of thickness ratio 1/10 at My = 0.8 are 
given in Table 3. Also included in Table 3 are the re- 
sults computed by Perl by the curvature method.’ 
Under the direction of the senior author, calculations 
have been made for compressible flow past bodies of 
As a fourth example, 
The complete 


The 


revolution by de los Santos. !'!” 
the case of flow past a sphere is taken. 
‘work of de los Santos will be published elsewhere. 
velocity potential ¢ may be assumed to be 


p _ o: + o- 
o; = [r + (1/2r?)| U cos 0 


l l l l 
Qe = Au ( a ee ) cos 0 + Aa ( eS :) cos 6 + 
2r? ty4 or* or 
l l 
Aj; op? 7 tr! cos* @ + Ao x 
l . l l 
— —]cos'§@+ As _— cos 6 + 
5r5 6r® Sr> 


l l ; 
Ajs (-. - =) cos’ @ (24) 


With ¢ assumed as in Eq. (24), the variational integral 


( 
373 


can be shown to be 
SSS (dn? — (0¢/dr)? — (Og/r00)2]~ r? sin 8X 
dr dd dw + 4yxrU (qn? — U2)” °°? [(Ay/3) + 
(Ai3/5) + (Aws/7)]/(y — 1) 


(25) 


where 7, 6, and w are spherical coordinates. Using y = 





FLOW 


EQUATIONS 347 











Y 
(Le 
3 
CIRCULAR CYLINDER 
Y 
Po - r bal Oe 
il 
4 
4 





/ 
/ 
/ o 
! P i 
: 
\ \ 
\ ~> 
\ 
\ 
\ 


_ 


“ 
~ 








~ =- 


—_——|—— 


ELLIPSE 


Y 


EE RS x 


t=.0513 








KAPLAN'S BUMP 


Fic. 2. 


TABLE 3 
Velocity Distribution g/U on the Surface of an Elliptic Cylinder 
of Thickness Ratio 1/10 at Mp = 0.8 


/ q Linearized Nonlinear Perl's 
U Solution Solution Solution 
0 0 0 0 
5 0.7174 0.6982 
10 0.9499 0.9256 
15 1.0270 1.0026 
20 1.0625 1.0399 
30 1.1006 1. O845 1.145 
40 1.1274 1.1199 1.1625 
50 1.1510 1.1526 1.1754 
60 1.1716 1.1819 1. 1827 
70 1.1878 1.2052 1. 1883 
80 1.1982 1.2203 1.1903 
90 1.2018 1.2255 1.1913 
TABLE 4 


Velocity Distribution g/U on the Surface of a Sphere at My = 0.50 


o/ q Linearized Nonlinear 
l Solution Solution 
0 0 0 
10 0.2473 0.2473 
20 0.4887 0.4881 
30 0.7206 0.7187 
40 0.9408 0.9376 
50 1.146 1.143 
60 1.331 1.328 
70 1.481 1.479 
80 1. 580 1.580 
90 1.615 1.616 








348 JOURNAL OF THE AERONAUTICAL SCIENCES 


2, the velocity distributions on the sphere at My = 
0.5 are given in Table 4. 

Comparing the results in these tables, we see that 
the linearized solution gives good approximation to the 
nonlinear solution up to the Mach Number close to the 
limiting value. Thus, for a given profile, if only the 
velocity and pressure distributions at various Mach 
Numbers are to be found, the linearized variational 
method offers a simple method for yielding these re- 
sults. The saving in labor in carrying out a linearized 
problem compared to that for a nonlinear problem is 
tremendous. For example, in the case of a Kaplan 
bump, by taking six parameters, some 540 terms appear 
in the final equations, but in the linearized case only 
42 terms are involved. The amount of saving in labor 
becomes obvious from a comparison in the numbers of 
terms. 
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Electrical Analog Solution for Centrifugally 


Tuned Pendulum 


Absorber System 


JAMES E. ANDERSON* ann WALTER W. SOROKAt 
Douglas Aircraft Company, Inc., and University of California at Berkeley 


SUMMARY 

The response of a two-mass free-wheeling dynamical system to 
sinusoidal excitation as determined experimentally by means of a 
simple electrical analog corresponds closely over a wide frequency 
range with theoretically predicted amplitudes for the mechanical 
system. The analog was extended to include a hypothetical 
engine-propeller system equipped with a centrifugally tuned 
pendulum absorber and subjected to several orders of engine 
excitation. The responses experimentally determined by this 
means agreed satisfactorily with theoretical calculations for both 
cases, with and without the absorber, for all harmonic orders 


present. 


INTRODUCTION 


Fe THE TIME OF ITS INTRODUCTION in 1936 by 
Taylor,' the centrifugally tuned dynamic absorber 
has been the subject of considerable theoretical analy- 
sis? and experimental work.* This, of course, is a direct 
result of the proved effectiveness of the device in sup- 
pressing troublesome orders of torsional vibration ex- 
citation in radial aircraft engines. This type of ab- 
sorber has also found application in the field of aircraft 
propeller vibration control.* 


In the application of the tuned pendulum to drive 
systems, it is necessary to consider its effect not only 
on the exciting order being suppressed but also on the 
other harmonic orders present in the drive system. 
The possibility of inadvertent magnification of the 
response at other orders of engine speed must be guarded 
against. This means that the dynamical characteris- 
tics of the entire drive system should be studied in 
detail as the pendulum tuning, mass ratio, and loca- 
tion are varied in order to arrive at optimum results. 
Moore® and Harker® have analyzed some of these prob- 
lems for simplified aircraft engine-propeller drive sys- 
tems. 


To design one or more absorbers and to specify their 
locations for optimum effects in complex drive systems 
is likely to become an analytical task of herculean 
The necessary computations may well be 
Of course, 


proportions. 
extremely laborious and time-consuming. 
direct experimentation may be resorted to, but the con- 
struction of the test samples, setting up of instrumen- 
tation, making provision for the variation of param- 
eters, and analyzing the records make this proce- 
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dure hardly practicable except for proving the final 
design. 

A promising experimental approach that may be 
conveniently and inexpensively undertaken in small 
laboratories, which also lends itself easily to wide varia- 
tion of parameters and to simple measurements, is 
that in which the electrical simulator,’ or analog com- 
puter, is used to model the mechanical system. Com- 
plex dynamical systems are as easily represented by 
this means as are the simple systems, the only extra 
requirement being the addition of more similar electri- 
cal elements. Many dynamical problems in vibra- 
tion and impact have been solved by such comput- 
ers.*~' Although large general-purpose analog com- 
puters'' '? have been set up in various laboratories 
(notably Westinghouse Electric Corporation and the 
California Institute of Technology), such machines are 
not readily available to the usual industrial engineering 
departments. It is entirely practical to set up small 
special-purpose computers in such engineering organi- 
zations for the solution of particular problems confront- 
ing them. 


THE MECHANICAL SYSTEM 


It is the purpose of this paper to show how a simple 
electrical model may be applied to the solution of vibra- 
tion problems involving centrifugally tuned dynamic 
absorbers and to present data for a hypothetical case 
that has been solved electrically and verified theoreti- 
cally. Since the analogy may be applied without sub- 
stantial increase in difficulty to problems of any prac- 
tical complexity, a simple case was chosen to illustrate 
the method. Fig. 1 is a sketch representing the sys- 
tem whose solution is presented herein. It might be 
considered as a first dynamical approximation to a 
radial engine-propeller system equipped with a pendu- 
lum vibration absorber. 


Vibration-exciting torques due to firing impulses and 
due to inertia effects in the engine reciprocating parts 
operate at frequencies that are directly proportional 
to engine speed. Consequently, the various excitation 
frequencies may be expressed as multiples or submulti- 
ples of engine speed. Taylor' shows that the gas- 
pressure exciting torques generated in a single-cylinder 
four-stroke-cycle gasoline engine consist of an infinite 
number of half-integral and integral multiples of engine 
speed, called harmonic orders (or engine orders or ex- 


349 











350 JOURNAL OF THE ABRONAUTICAL SCIENCES 























Fic. 1. 
cal analog (bottom). 


citing orders). Fortunately, the exciting torque am- 
plitude drops off rapidly with increasing engine order, 
so that only the first few orders need be considered in 
any practical analysis. In multicylinder engines, se- 
quential firing tends to reduce or to eliminate many of 
the harmonic frequencies, usually leaving as the most 
important order the one equal to the firing frequency. 
Wilson'* tabulates typical amplitudes of the harmonic 
orders for single-row, seven- and nine-cylinder radial 
engines. Predominant exciting frequencies are the 
three and one-half order for the former and four and 
one-half order for the latter. However, appreciable 
first- and second-order harmonics exist in both cases, 
as well as an important seventh-order harmonic for the 
first engine and a ninth-order harmonic for the second 
engine. 

Bentley and Taylor": present data for the gas-pres- 
sure torque harmonics of two-row radial engines. By 
arranging crank throws at 180°, the seventh- and ninth- 
order harmonics are substantially eliminated for the 
14- and 18-cylinder engines, respectively. Varying 
the relative master cylinder positions through all pos- 
sible settings permits wide adjustments in first- and 
second-order harmonic amplitudes. Misfiring intro- 
duces a considerable half-order harmonic component, 
as well as affecting appreciably many of the other har- 
monics. 

It is evident, then, that the designer of a complex 
drive system (which necessarily includes the engine) is 
faced with a difficult problem of determining the critical 
frequencies of the system, the expected responses in 
resonance with the various exciting orders likely to be 
present, and of providing a means for reducing the 
critical responses to tolerable limits. Generally, seven- 
and nine-cylinder, single-row radial engines (and the 
corresponding 14- and 18-cylinder, double-row en- 
gines) are equipped with three and one-half and four 
and one-half order pendulum dampers, respectively, 
since these are likely to be the most serious exciting 


Representative dynamical system (top) and its mechani- 


JUNE, 1950 


orders. However, serious vibration difficulties may 
exist at some other harmonic orders as well and may 
necessitate centrifugally tuned absorbers to reduce ey. 
citation from these sources. 
absorber‘ provides such an example. 

In the hypothetical case presented in this paper, the 


The propeller vibration 


important harmonic orders are considered to be the 
third, fourth, fifth, and sixth, with the fourth order 
predominant. Consequently, a fourth-order pendulum 
absorber is provided, and its influence on the system 
response is determined. A theoretical analysis of the 
system with and without the absorber has been made, 
tests were run on an electrical model, and the two sets 
of results are compared. 

Without the absorber, the equations of motion of the 


system of Fig. 1 are 


1,6; + 0,6, + (A; — hh) = T\ (1) 


156» + bb. + k(05 —_ 6) = 0 f 
where 
J, = equivalent engine inertia 
I, = equivalent propeller inertia 
k = equivalent stiffness of propeller shaft 
b, = viscous damping coefficient associated with the 
engine 
b, = viscous damping coefficient associated with the 
propeller 
7 = one of the harmonic orders of the exciting 
torque 
Let 7 = Tye’. Whereupon the Eqs. (1) have the 


steady-state solution 


A; = Ad", Oo = Ave™ 
where 
A, = To(k —_ Toa)” + lw») [(Rk ag Tw? + 
twh,)(k —Isw* + iwh.) — k*I 


Ay = Tok/[(k — Tw? + twh,)(k — Lew? + twhe) — k’| 


The absolute amplitudes of response are, therefore, 


4) = TH — oh — FO = #Id/0- Ah 4 
wb»*)| ? + [w, + k*whe (k 2 “ie oy w%by?)) 2)" (9) 


— Tok \[(k si Tw") (k = Tou") me whi bo — k?|? + 

w*[bo(k — Iw?) + d(k — Tow?)|?} (3) 

The following numerical values were assigned to the 
system properties: 


I, = 10 in.Ibs. sec.” a = S48 in.Ibs. sec. 
I, = 50 in.Ibs. sec.? b, = 2,390 in.Ibs. sec. 
k = 4X 10° in_Ibs. 7) = 6,000 in.Ibs. (third 
r.p.m. = 1,000-2,500 order) 


15,000 in. Ibs. 
(fourth order) 
= 5,000 in.Ibs. (fifth 
order) 
= 7,000 in.Ibs. (sixth 
order) 


Il 
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Cc = capacity of the condenser, farads 


TNE R,, Ro = resistances, ohms 


rw a : 
a oe @, Q@ = electrical charge, coulombs 
In going from the mechanical elements to the analo- 
a gous electrical elements, the following relations were 


a a used : 5 














L=al, q= 06, t. = ct (5) 

, genoa where a, ), and ¢ are the conversion factors and ¢, and 

L R L R. t, are the times (seconds) for corresponding events to 
« . fe occur, respectively, in the electrical and mechanical 


wc. 2. Electrical analog of two-mass free-wheeling system. se . : 
a . ie systems. [Ep. Note: No Eq. (6).] 
In terms of the fundamental conversion factors a, 3, 


A detailed calculation of the response of the engine sie ; ; 
and c, the remaining derived conversion factors are 


inertia J; to the fourth-order exciting torque was made 
over a wide range of engine speeds, including the entire E = (ab/c?)T 1 = (b/c)w } 
normal operating range of the engine. The theoretical R, = (a/c)b, L/G (a/c?)k ? (7) 
response iS represented by the full-line curve of Fig. 3. R, = (a/c)bo \ 

Practical magnitudes for the electrical circuit ele- 
ments were obtained by taking the following values for 


THE ELECTRICAL ANALOG 


The system without the absorber was simulated by — the fundamental conversion factors: 

the electrical circuit of Fig. 2, making use of the mass- ee . r 
_ ‘ . ; a = 2.45 X 10~‘4 henrys per in.Ibs. sec.” 
inductance analogy. In Fig. 2, an evaluation of the ang 1 ‘ . 

‘ : . b = 5.77 X 10~* coulombs per rad. (8) 
sum of the voltage drops around each loop in accord- a ee : \ 

~ : : c = 10~* (dimensionless) 
ance with Kirchhoff’s law yields the equations 

The magnitudes of the elements in Fig. 2 are, there- 


Ligs + Rigi + (1/ Om — @) = El (4) oe 
Loge si Roge + (1 C) (qe —- fA) = 0 f j 
L, = 2.45 millihenrys, L. = 12.25 millihenrys 
ee C = 0.102 microfarad, EH = 0.212 volt (peak) pp (9) 
L;, L, = self-inductances of the coils, henrys R, = 20.8 ohms, R: = 58.6 ohms \ 
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Fic. 3. Frequency-response curves for two-mass free-wheeling system 
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Fic. 4. Method for exciting analog circuit. 
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The resistance values R,; and R, given above are nom- 
Actually, these magnitudes varied with 
frequency. Over the test range, the values ranged 
between 19.7 and 21.5 ohms for R; and between 57 and 
60 ohms for Ro. 

The electrical system is 100 times as fast as the 
mechanical system by virtue of the choice of c. Thus, 
the natural frequency of the electrical analog is 100 
times that of the mechanical system. Similarly, the 
forcing frequencies impressed on the electrical system 
must be 100 times the corresponding forcing frequen- 
cies in the mechanical system. 

An audio oscillator having a peak output of approxi- 
mately 30 volts root-mean-square provided the input 


inal values. 


Circuit to determine effective a.c. resistance of analog 


Electrical model of tuned pendulum engine drive system. 
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signal, #. The oscillator output was impressed acrogg 
two resistances in series (Fig. 4), the resistance R,’ 
being small compared with both R,’ and the impedance 
of the analog circuit proper. The exciting signal 
was picked off the terminals of R,’. 
ment permits accurate control of the input FE and also 


Such an arrange. 


reduces to negligible proportions load variations jm. 
posed on the oscillator by wide variations of impedance 
in the analog circuit with changing frequency. 

The inductance coils employed (shown in the photo. 
graph of Fig. 6) are commercially available radio-fre. 
quency choke coils wound so as to minimize distrib. 
uted capacitance. They had some resistance; hoy. 
ever, additional resistance was provided in series with 
the coils in order to supply the specified amount of 
damping. The total d.c. resistance in the circuit was 
considerably less (R; = 9.3 ohms, R2 = 20.8 ohms) than 
the a.c. resistance in operation. Since the a.c. resist- 
ance was the important magnitude, each coil and its 
associated resistor were checked on a separate test for 
a.c. resistance over the frequency range of the test. 
The checking setup is shown in Fig. 5. A variable 
condenser, C’, was used to bring the circuit into reso- 
nance at a series of frequencies over the desired test 
range. At forcing frequencies other than the reso- 
nance frequency, the current in the loop is given by 


T = E,/[Rey? + (Xi — Xc)"I' 
At resonance X,; = X¢, whereupon 
Dies. = By Rosy (10) 


If one were to measure /,,,., the above equation may 
be used to determine R,,, at the given frequency. In 
the present test, it was simpler to measure £, across the 
condenser C’ and to make use of the resonance relation 
X, = Xc in order to write 


FE, = Xclree. = Xt res. = QafLT res. = 2afL(Es/Re) 


whence, 
Ruy = 2afL (ke Ff) (11) 


Variations in capacitance and self-inductance over 
the frequency range of the test were considered negli- 


gible. The decade condenser was rated by the manu- 


facturer as having less than 0.25 per cent variation in 
capacitance over this frequency range. Also, its dis- 
sipation factor was extremely low over this range 
about 0.0005. Since the inductance coils are air-core 
coils, their inductances are substantially constant over 
the test range. The distributed capacitance of the 
coils and possible self-inductance effects in the resis- 
tors R; and Ry were neglected. The good agreement 
between the theoretical curve of Fig. 3 and the exper- 
mental points spotted thereon indicates that these ef- 
fects were negligible. 

The next step was to add to the analog in Fig. 2 4 
loop representing the tuned pendulum absorber. In 
order to establish appropriate values for these circuit 
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CENTRIFUGALLY TUNED 


it was necessary to go back to the theory of 


elements, 
the pendulum absorber for values of equivalent stiffness 


The theory has been excellently de- 
veloped in literature.” > ° For present purposes it 
is sufficient to indicate that the equivalent stiffness 
k, of the point-mass pendulum, m, on its suspension in 
a centrifugal field of angular velocity 0 is 


and inertia. 


k, = m(R/r)(R + r)?Q? (12) 
where 
R = radius from shaft center to pendulum attach- 
ment point 
r = length of pendulum 
m = mass of pendulum 


The inertia of the pendulum is 


I, = m(R + r)? 


Pp 
Its natural frequency about its pivot is 


Wyn = (k,/I,)”* = Q(R/r)” (13) 


Eq. (13) points out the characteristic property of the 
pendulum in a centrifugal field—namely, that its nat- 
ural frequency is proportional to the speed of rotation, 
Q. The order of tuning of the pendulum is determined 
by the ratio (R/r)”*. 

Viscous damping in the pendulum absorber is ex- 
tremely low. This is in contrast with the fixed-fre- 
quency tuned absorber, where a certain optimum 
amount of damping is necessary in order to level out 
response peaks over the operating speed range. Conse- 
quently, the pendulum absorber has been considered 
frictionless in deriving Eq. (13). 

The commercial radio-frequency coils available in 
the laboratory were considered to have too much re- 
sistance to represent the absorber mass. A _ special 
Brooks-type coil'® was wound of relatively heavy wire. 
In conformity with the chosen conversion factors, the 
absorber coil inductance was adjusted to a value of 
0.443 millihenry, corresponding to a mass, m = 0.025 
lb.-sec.? per in. at an equivalent distance of 8.5 in. 
from the shaft center. Since a fourth-order damper was 
to be investigated, the equivalent pendulum length was 
set at 0.5 in., with the pivot point 8 in. from the shaft 


center. Thus, 


Wp,/Q = (R/r)”? = (16)? = 4 


The capacitance, C,, corresponding to the equivalent 
flexibility of the absorber connection is determined 
from 


C, = (b?/a)(1/k,) (14) 

Since k, is a function of the speed of rotation, the 
value of the capacitance C, must be adjusted to its 
proper value for each excitation frequency. Conse- 
quently, a decade condenser was used for this pur- 
pose. For each setting of the condenser, the oscillator 


Was adjusted to provide the appropriate excitation fre- 
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Electrical analog of dynamical system with pendulum 
absorber. 


Fic. 7. 


quency and amplitude in accordance with the particu- 
lar order of vibration under investigation. A photo- 
graph of the actual electrical model is shown in Fig. 6. 
The corresponding schematic circuit is shown in Fig. 7. 

The measured responses of the engine mass /, (con- 
verted back to mechanical units) are shown plotted in 
Fig. 8 for the various exciting orders in the absence of 
the vibration absorber and in Fig. 9 for the same condi- 
tions with the absorber present. A comparison of the 
two figures shows strikingly how the absorber affects 
the amplitudes of torsional vibration under the various 
orders of excitation present in the system. The fourth 
order almost disappears. The other orders are 
greatly affected in amplitude for this particular case, 
but there is a considerable shift in the speed at which 


not 


maximum responses occur. 

As a direct check on the accuracy of the results, the 
theoretical frequencies and amplitudes of maximum 
response were calculated for the two cases, with and 
without the absorber. This calculation is simplified 
by virtue of the development brought out clearly by 
Stieglitz’—that the effect of the absorber in each ex- 
citing order is equivalent to that of an increase in in- 
ertia of the disc to which it is attached, the added 
amount being dependent on the order of excitation but 
not on the actual speed of operation. This increase is 
given by 
(15) 


AI = mR(R + 1r)?/(R — rn?) 


where » = w/Q2 = harmonic order of excitation. This 
is an important property of the tuned pendulum ab- 
sorber, in that the introduction of the absorber does 
not bring in with it a spring and, therefore, the addi- 
tional resonance that is characteristic of fixed-frequency 
tuned dynamic dampers. 

Thus, in calculating the response of the system, it is 
appropriate to replace the inertia /; with an equivalent 
inertia given by 


| Oe = I; + Al (16) 


which is to be evaluated for each order of excitation. 
Eq. 2, previously derived for the two-mass system, may 
now be used with /;,, to evaluate the theoretical response 
of the system with absorber. Since the resonant fre- 
quencies are not critically dependent on damping, they 
may be calculated without including damping by means 
of the following expression :'7 
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Engine response curves without absorber. 


LEGEND: 
© CALCULATED Points 


-eo—e EXPERIMENTAL CURVES 


3RD ORDER 


28 30 32 34 36 38 ” 


ENGINE SPEED (res) 
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f, = 30 [RU Lieg + Is) 
; 7 N Dregl 


cycles per min. 
The critical speeds of the various orders will be 


Nerit = Se n 


Table | correlates values of the calculated and meas- 
ured frequencies and amplitudes of the resonance peaks 


for the various orders. The experimental and theoreti- 


Engine response curves with fourth-order absorber. 


cal results agree to about | per cent in the frequency 
determinations. The maximum deviation betweet 
theory and experiment in resonant amplitude deter 
minations is 7 per cent. Since resonant amplitudes are 
critically dependent on the amount of damping pres 
ent, the greater error in peak amplitude measurements 
is due to the greater uncertainty in assessing the resist 
ance of the electrical circuit. With more precise ele 


ments it should be possible to reduce this discrepancy 
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CENTRIFUGALLY TUNED 
TABLE | 


Resonance Frequencies and Amplitudes of Engine Mass 
Without Aksorber 


Com- Experi- 
puted mental Com- Experi- 
Torque Fre- Fre- puted mental 
Ampli- quency, quency, Ampli- Ampli- 
Engine tude, Cycles Cycles tude, tude, 
Order In.Lbs per Sec. per Sec. Deg. Deg. 
8rd 6,000 36.7 36.8 0.527 0.515 
4th 15,000 27.5 27 .6 1.318 1.267 
5th 5,000 22.0 22.1 0.439 0.430 
6th 7,000 18.3 18.4 0.615 0.590 
Resonance Frequencies and Amplitudes of Engine Mass with 
Fourth-Order Absorber 
Com- Experi- 
puted mental Com- Experi 
Fre- Fre- puted mental 
quency, quency, Ampli- Ampli 
Engine Cycles Cycles tude, tude, 
Order per Sec. per Sec Deg. Deg. 
3rd 32.0 32.1 0.579 0.550 
5th 26.0 26.0 0.381 0.355 
6th 19.6 . 19.8 0.563 0.530 
System Properties: 
I, = 10 in.lbs. sec.2 R = 8.0 in. 
I, = 50 in.lbs. sec.2 r = 0.5 in. 
k = 4X 10%in.lbs. m = 0.025 lb. sec.” per in 
considerably. However, it is difficult to assess damp- 


ing in engine installations, so that it is doubtful whether 
the damping coefficients for the mechanical system 
would be known within 20 per cent for an actual in- 
stallation. Consequently, the experimental amplitude 
determinations are considered to be entirely adequate 
for the exploration of absorber designs and locations. 
The final design should be checked in full-scale tests. 


CONCLUSION 


The theoretical and experimental work performed in 
connection with the design of a hypothetical centrifu- 
gally tuned pendulum absorber shows that an electrical 
model may be conveniently and economically used 
in the design and analysis of complex systems requiring 
With reasonable 
care accurate results may be obtained. The long and 
involved calculations ordinarily necessary in a numeri- 
cal analysis may be replaced with measurements on the 


the application of such absorbers. 


electrical model. 
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The Shearing Rigidity of Buckled Sheet 


Panels 


BRUNO A. BOLEY* 
Goodyear Aircraft Corporation 


SUMMARY 


A formula giving the effective shearing rigidity of either flat 
or curved panels in the buckled state is developed from approxi- 
mate theoretical considerations making use of the concept of 
effective width. Cases of panels buckled by compression, shear, 
or combined compression and shear are considered. The formula 
contains an experimentally determined exponent and gives values 
both for an effective shear modulus to be used in determinations 
of stress distributions and for a reduced effective shear modulus 
to be used in instability calculations. An experimental check is 
established. 


INTRODUCTION 


ymin OF THE Shearing rigidity of panels of 
thin sheet is often required in the analysis of 
stressed-skin constructions, both for the determination of 
stress distribution (as in shear lag analyses) and in in- 
stability calculations for sheet-stringer combinations 
or monocoque cylinders. In the present paper, an ap- 
proximate formula is developed which gives the shear- 
ing rigidity of sheet panels, as represented by an ef- 
fective shear modulus, for simple, as well as combined, 
loadings and for both flat and curved panels. 

If a panel is not in a buckled state, its rigidity is rep- 
resented by the shear modulus Gp of its material. The 
relation between the shear deformation and the stress 
that causes it is then Hooke’s law: 


t = Goy (1) 


where 7 is the shear stress and y is the shear strain. 

If the panel is in a buckled state, however, the de- 
formations resulting from the application of a given 
shear stress will depend on the complex state of dis- 
tortion of the panel. Eq. (1) is therefore no longer 
valid, but an equation of the same form as Eq. (1) will 
still hold, provided that an effective shear modulus of 
Gi. be used in place of Go. However, since the shear 
stress and strain vary from point to point in a buckled 
panel, it is necessary to make use of an average value 
of the shear stress 7,,,. and of the shear strain ya), 
Then 


Tare. = Goss. Vase. (2) 


and G,y,, = Go if the panel is not buckled. For buckled 
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panels G,yy, is lower than Gp because of the increased 
movability afforded by the waves or folds. Thus, this 
reduction in the shear modulus is not due to an inelastic 
behavior of the material, and, in fact, any element of 
small dimensions in the panel will, in general, follow 
Hooke’s law even after buckling. The rigidity of the 
panel is nevertheless reduced because part of the shear. 
ing deformations are taken up, within the panel, by 
straightening the folds and reducing the slack that the 
buckling waves introduced. This action requires 
only bending and no stretching of the sheet and is met 
with little resistance because of the small thickness oj 
the material. It may then be concluded that the shear. 
ing rigidity of a buckled panel depends on both the 
amount of waviness present (which is a function of the 
amount of load) and the shape of the waves (which isa 
function of the type of loading). 

The effective modulus defined by Eq. (2) may be 
used in calculations of stress distribution. In problems 
of instability, however, it is desired to know how much 
resistance the structure will offer against distortions 
additional to those existing immediately before the in- 
stant of buckling. If the existing distortions are rep- 
resented by the shear strain 7Yq,- (corresponding to the 
stress T,,-.), then the additional distortions may be rep- 
resented by dyq,.. (corresponding to an additional stress 
dtare.). The relation between the additional shear 
strain and stress will again have the same form as Eq. 
(1), provided that a reduced effective shear modulus 


Ge sy. rea. be used in place of Gp. Then, 


Or ese. = Gory red dV ave (3) 


Similarly to the modulus G,,,, the modulus Gyyy._ rea will 
also depend on the amount and type of loading on the 
panel. It was first used in reference 1 and was shown 
there to be, in general, lower than G,,,. 

Values of the effective shear modulus have been ob- 
tained by many authors, some of whom are listed here. 
Results obtained with flat panels under shear are pre- 
sented graphically in reference 2, while a theoretical in- 
vestigation of flat panels under combined loading led to 
the curves of reference 3. The torsional stiffness of 
circular cylinders was investigated experimentally i 
reference 4 for circular cylinders under no compressive 
loads and in reference 5 for cylinders under various 
amounts of compression. In the latter report, e* 
perimental curves were presented, as well as an em- 
pirical formula in which the effective modulus de- 
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REGIDSDiITyY oF 


SHEARING 


pends on the angle subtended by the panel and on the 
ratio of the applied strain to the buckling strain of the 
panel. These two parameters were established by 
physical reasoning in reference 6. A theoretical anal- 
ysis of flat panels under shear above the buckling load 
was carried out in reference 7 and led to values for an 
effective width of panel corresponding to the effective 
shear modulus. Values for circular cylindrical panels 
under combined compression and shear were needed for 
investigations such as those of references 1, 8, or 9, 
but they could not be found in the literature. 

This brief review of past investigations shows that 
there is a need for determining more values of the shear 


BUCKLED 
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physical action of the buckled panels have been kept. 
Furthermore, the usefulness of a more accurate analysis 
is impaired by the many uncertain factors encountered 
in actual problems of buckling of thin panels (such as 
end-fixity conditions or initial distortions). The formu- 
las contain an empirical exponent, but, since a single 
constant value for this exponent gives consistent re- 
sults for all cases of loading it may be said that the 
formulas are in good agreement with experiment. 

The formulas for the various loading conditions are 
collected in Table 1. 
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modulus for some of the loading conditions, as well as a = side length of panel 
for devising a convenient method for obtaining the b ne width of panel 
for all loading conditions without referring to be, be ° Geen oe 
modulus for a . 8 é 8 E = Young’s modulus 
a large number of graphs in several separate reports. Geff. ws eieten deus unetien 
The approximate formulas developed in this paper on Geff. red. = reduced effective shear modulus 
the basis of von K4rman’s effective width concept are Go = shear modulus of panel material 
believed to answer these needs and to be sufficiently L ms icici 
accurate for most purposes in design and stress analysis. 4 i a 
While the development of the formulas is obviously not t é Ghat tee 
rigorous, it is believed that the essential features of the y = shear strain 
TABLE 1 
FORMULAS FOR THE SHEARING RIGIDITY OF SHEET PANELS 
Loading Type of | Effective Shear Modulus | Reduced Effective Shear Modulus | When 
anel Gott Gers. red. valid 
G Gq Ys G, -red. G. 
Shear Flat te) et = 2 a 
Go Tave G, ° 
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Yave. = average shear strain ie 
v = Poisson’s ratio 
o = compressive stress UH L Li} 
cr. = critical compressive stress under combined loads 
ocu. = compressive buckling stress of curved panel [given 
in Eq. (21)] 
Tu. 0 = compressive buckling stress of unreinforced circu- 
lar cylinder [given in Eq. (22) ] a ' 
of = compressive buckling stress of flat panel [given in Y = — 
Eq. (9) | 
T = shear stress Wi | | | iil 
Tave. = average shear stress 
Tor. = critical shear stress under combined loads 
Tou. = buckling shear stress of curved panel [given in Oo {I Oo 
Eq. (21)] 
Tou. 0 = buckling shear stress of unreinforced circular cyl- was b 1 L be 
inder [given in Eq. (20)] Be el 
Tf = buckling shear stress of flat panel [given in Eq. 2 2 
(10)] 


THE FORMULA FOR THE EFFECTIVE SHEAR 
MopuULus 


A flat panel of width 6 subjected to uniform edge- 
wise compression will be considered first. If the panel 
is in a buckled state and its edges are reinforced by 
straight bars, the regions near the edges of the panel 
will be more rigid than the unsupported middle por- 
tion. The load-carrying capacity of the panel may 
therefore be assumed, as is often done, to be concen- 
trated in two strips of total width b,, along the edges of 
the panel, as shown in Fig. la. 

Now let the panel be subjected to a shearing action 
as shown in Fig. lb. The effective width strips may 
be assumed to follow Eq. (1). From Fig. 1b, 


Y = Yave. (4) 
The total shear force F applied to the panel is 
F = chg = Tas. dt (5) 


from which the average stress may be calculated. Sub- 
stitution into Eq. (1) gives 


Tou, = Go (b, b) Yave (6) 
and comparison with Eq. (2) yields 
Gory, ‘Go = b, b (7) 


From this formula the effective modulus can be calcu- 
lated if the effective width b, is known. The deter- 
mination of the effective width will be considered later, 
and the above approach will first be applied to panels 
buckled by loadings other than pure compression. 
Consider the combined loading of Fig. 2a. The buckles 
will extend in the general direction of the tensile (or 
smaller compressive) principal stress, while the effective 
width strips will extend in a direction approximately 
perpendicular to that of the buckles. This situation 
is indicated in Fig. 2b. Now, again, it may be as- 
sumed that only the effective width strips have load- 
carrying capacity, so that it is necessary to consider 
only the effective panel of Fig. 2c, which can, as before, 
be assumed to follow Eq. (1). A development entirely 
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Fic. 1. Panel buckled by compressive stresses 
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Fic. 2. Buckled panel by a general type of loading 


analogous to that used with panels under compression 
will again lead to Eq. (7). Moreover, this development 
can be easily extended to curved panels. Eq. (7) 1s 
therefore suggested for the calculation of the effective 
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Fic. 3. Deflections in the effective width strips 


modulus of thin panels, the effective width b, being dis- 
cussed in the next section. 

It should be pointed out that the shearing rigidity 
obtained from this approximate theory will be some- 
what smaller than the actual one, since, in reality, the 
middle portion of the panel, the rigidity of which is now 
neglected, also makes some contribution to the total 
stiffness of the panel. However, the use of an em- 
pirical exponent may be considered to correct this 


error. 


THE EFFECTIVE WIDTH 


An approximate method for calculating the effective 
width of a flat panel under compression has been de- 
veloped in reference 10. This method will be used now 
to obtain the effective width of panels under other 
loadings. 

Let the curve ABC in Fig. 2a represent the deflected 
shape at the instant of buckling of line DE, which is in 
the direction of the buckles. Point B is the point of 
maximum deflection, so that the plane tangent to the 
deflected plate at that point is parallel to the plane of 
the unloaded plate. It may be assumed in first ap- 
proximation that above the buckling load the portion 
of line DE which is contained within the effective width 
strips is deflected according to curves AB’ and B’”A 
(Fig. 3) of the same shape as curve ABC but, of course, 
shrunk into a smaller length. The portion of the panel 
outside of the effective width strips may be assumed to 
be straight. Then, the effective panel (Fig. 2c) will be 
under exactly the same deflected shape (of course, to an 
appropriate scale) which existed in the complete panel 
at buckling. The stresses o,,. and 7,,, carried by the 
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complete panel at buckling under this deflected shape, 
are given by the equation 


(Ger./ 6s) + (Ter./ts)? = 1 (8) 


according to reference 11. The buckling stress o, of a 
flat panel under compression alone is shown in refer- 
ence 12 to be 


oy = [ki w? E/12(1 — v?*)] (t/b)? (9) 


while the buckling stress 7, of a flat panel under shear 
alone is, from the same reference, 


ty = [ke x? E/12(1 — v?)] (t/b)? (10) 


In these equations, k; and k2 are nondimensional fac- 
tors dependent on the ratio (a/b) of the sides of the 
panel and on the end fixity of the edges and may be ob- 
tained from reference 12, pages 330 and 360, respec- 
tively. 

The stresses o and 7 carried by the effective portion 
of the panel can be obtained, following the previous 
discussion, from the above formulas, provided that the 
effective width b, is used in place of the total width bd. 


Thus, 
o (j) 4 E (i) | ue (11) 
ar \b tT, \b 


The effective shear modulus is then given by 


@ (Sen) + (=) (4) —— (12) 
os Go Ty, Go 


For the special case of pure compression, this reduces 
to the well-known von K4arman’s formula 


b./b = Va;/o (13) 
It is known that this formula gives values that are 
smaller than the corresponding experimental results 
when applied to the calculation of the width of panels 
effective in carrying compressive stresses. In this 
case, Marguerre’s formula'* 

b,/b = (a,;/0)* (14) 


is commonly used in place of Eq. (13). In the present 
analysis, it will be found that a similar change in the 
numerical value of the exponent will be necessary to 
bring the results into agreement with experiment. 


Therefore, Eq. (11) will be written as 


1/n 2 2/n 
o ("" n (7) (** ~ rs 
a, \b tr) \b 


and the value of the exponent will be discussed later, 
together with the empirical check. 

It is more convenient to obtain the effective width 
(and the effective shear modulus) in terms of the aver- 
age shear stress T,,- rather than the stress r carried by 
the effective width. This may be done with the aid of 
Eq. (5). The effective width is then the solution of the 


equation 
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1/n 2 2(1—n)/n 
.. (;*) 4. (=) (;') at (16) 
a, \b Ty b 


which can be solved explicitly for certain values of n. 

The effective width of a curved panel may be ob- 
tained by adapting a suggestion of reference 14. Let 
T-. be the shear stress that the panel can sustain because 
of its curvature; then, the total shear force F may be 
written as 


F = Tare. bt = Tey, Ot + (Tare. — Tex.) Ot (17) 


in which the last term is the portion of the force car- 
ried by the panel independently of curvature. This 
portion is carried by an effective width b,’, which may 
be calculated from Eq. (16), provided that the stress 
(Tave. — Teu.) be used in place of the stress Ta,-. appear- 
ing in that equation, and, similarly, for the stress o. 
Thus, 


(: ex, cen (%) re (= = ray n = ; 
os b uf b 


(18) 
Asa first approximation, it is possible to write 


b.’ = b (19) 


€ 


where b, is the total effective width. This means that 
the rigidity provided by curvature is considered not 
directly, but only inasmuch as it decreases the stress to 
be carried without the aid of curvature. This assump- 
tion leads to a shearing rigidity that is smaller than the 
actual one and is therefore on the safe side. The error 
introduced by it is, moreover, believed to be small. 

The stress 7,,. may be taken as the buckling stress 
of an unreinforced circular cylinder. If no compressive 
stresses are acting, 


‘oa. = 


(O.8S15E(r/L)"* (t/r)"*, if (L/r) < 6.5 (r/t)” 


= < 


(0.29 E(t/r)* (otherwise) (20) 


Teu. 0 
according to reference 15. If compressive stresses are 
present, the formula'® 


9 ‘ 
(Teu. Teu.o)” = B- (eu. Feu. 0) (21) 


may be used, where the buckling stress o,,. 9 of an un- 
reinforced circular cylinder under uniform axial com- 
pression is 


Cou. 0 = 0.3 Kit Yr) (22) 


The buckling stress o,,. is given by Eq. (22) if buckling 
is caused by compression alone. If both shearing and 
compressive stresses are acting, o,,. is given by Eq. (21) 
in terms of 7,,... Since the evaluation of the effective 
width from Eq. (18) requires the values of both o,,,. and 
Tey.» another relation between these two quantities is 
needed. It is reasonable to assume that the shearing 
and compressive stresses carried because of curvature 
are proportional to the total applied shearing and com- 
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Fic. 5. Effective shear modulus for flat panels buckled by 
compressive stresses. Suggested value for m = '/, 


pressive stresses, respectively. Then, 


T Teu. = o Feu. (23) 


A combination of Eq. (7), (18), and (19) gives the 
desired formula for the effective shear modulus as 


o—o ) G vies Ten. — Ta 1G" 
cu. esf. +4. ( ave ) ( u:) = | 
( Oo; ( Go ) Ty Go 


(24) 


COMPARISON WITH EXPERIMENT 


Values of the effective shear modulus given by the 
above formulas are plotted in Figs. 4—7 for the limiting 
cases ¢ = 0 andr = 0 for both flat and curved panels. 
Experimental results are also given as obtained from 
references 2 and 5. The experimental results shown 
in Fig. 6 were obtained with a cylinder described in 4 
thesis submitted at the Polytechnic Institute of Brook- 
lyn by Adolph J. Athen, whose contributions the au- 
thor gratefully acknowledges. The principal char- 
acteristics and dimensions of this cylinder were as fol- 
lows: 
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Radius, ry = 8 in. 

Sheet thickness, ¢ = 0.012 in., 24S-T Alclad 
Stringer section, */s by */s in., 24S-T 
Aluminum Alloy 

Ring section, */s by */s in., 24S-T 
Aluminum Alloy 

Number of stringers, 16 

Number of rings, 6 

Ring spacing, 4.5 in. 

Length of cylinder, 31.5 in. 


Equal and opposite torques were applied to each end of 
this cylinder through heavy steel end rings. The rela- 
tive angle of twist between the ends of the cylinder was 
measured by means of two Ames dial gages and plotted 
against the total applied torque 7. The average shear 
stress was obtained from the formula 


Tare. _ i (2A ,t) 


where A; is the area included by the median line of the 
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skin. The average shear strain y,,,. is related to the 
angle of twist ¢ by the equation 


Yave ” gr L 


where L is the length and r is the radius of the cylinder, 
so that the effective shear modulus is 


Geyy. _ Tare Yare = FB (2A jtre) (25) 


By multiplying numerator and denominator by the 
circumference S = 2zar of the cylinder and by rear- 
ranging terms, the well-known Bredt’s formula?’ 


Gey, = TLS/(4A 7te) (26) 


is obtained. 

The values of G,,y, given by the formulas previously 
derived are presented in Figs. 4-7, for values of the 
exponent m = (1/2), (1/3), and ('/4) [see Eqs. (7) 
and (15)]. It may be noted that the curves for n = 
'/, are in best agreement with the experimental results 
in Figs.4,5,and7. In Fig. 5 the shear moduli obtained 
on the basis of » = '/, are substantially lower than the 
empirical values. For comparison, curves are included 
in Fig. 5 which correspond to the theoretical results ob- 
tained in reference 3 under various edge conditions. 
They show that, in this case, the shear modulus varies 
greatly with edge conditions. As the edge conditions 
are difficult to determine in an actual case, the curve 
corresponding to nm = '/4, which is an average between 
the extreme values given, may be considered satisfac- 
tory. 

It is therefore suggested that the value n = '/, be 
used in the calculation of the effective shear modulus. 
It may be remarked that the value m = '/, gives curves 
that represent a good average of the experimental points, 
while n = '/; gives values that are always on the safe 
side. More experimental data are needed, however, 
especially in cases of combined loadings. The formu- 
las corresponding to n = '/,4 for various loadings and 
for both flat and curved panels are collected in Table 
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1. Values of the effective shear modulus for combined 
shear and compression, as given by the appropriate 


) 


formula of Table 1, are plotted in Fig. 8. 


THE REDUCED EFFECTIVE SHEAR MODULUS 


The reduced effective shear modulus G,,,. ;.4, was de- 
fined in Eq. (3). In this equation, the quantity 
(dT are.) Was defined as the change in shear stress due to a 
change in shear strain (d7Yqy-e.)._ Thus, in the calcula- 
tion of the reduced effective shear modulus, it is not 
necessary to consider any variation of 7,,,. with the 
compressive strain. The reason for this is that in in- 
stability problems Gy, ;.¢. is used in the evaluation of 
the work done by the shear stress, which is, of course, 


zero if only a compressive strain and no shear strain is 


present. Eq. (3) may then be written as 
a) Fass. d( Yave G, SS ) ‘ d G, ey. 
Gers. red. — = = Ger. + Yave. 
dar dar dave 
(27) 


If the effective modulus G,,,, is known in terms of 
stresses rather than strains, a similar formula may be 
obtained as follows: 


Gers we = = 


OVae Oren d( Tave Gory ) Tass 
G 
eff 92 
: ae (28) 
l (Tave./Gery.) (dGeyy./AT ave.) 
As was shown in reference 1, it follows from these equa- 
tions that 


(1) Gegy. rea. = Go if the panel is not in a wrinkled 
state. 
(2) Gery. rea. = Ger. if the average shearing strain 


(or stress) in the panel is zero immediately before 
buckling. 


Gey red 


According to the previous discussion it is suggested 
that n be taken equal to '/,. 
las for various loadings and for both flat and curved 
panels are collected in Table 1. 


The corresponding formu- 
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(3) Gegy. rea. < Ges. in all other cases, since, in gen 
eral, the effective modulus G,,, decreases with increas- 
ing strain (or stress) so that the terms dGyyy/d Yq... and 
IG eyy./d Tare. are Negative. 

With the aid of the formulas previously derived for 
the modulus G,,,., it is possible to eliminate the deriva- 
tives appearing in Eq. (27) or (28). It is necessary to 
consider only the case of curved panels under com- 
bined compression and shear, since all the other cases 
investigated are special cases of this. The modulus 
Gey. rea. May then be obtained by differentiating Eq. 
(24) with respect to 7,,,., solving the resulting expression 
for dGyzy./dTaye. and substituting the result into Eq. 
(28). If use is again made of Eq. (24), the following 
formula is obtained: 


© — Oey.\{Gery ales 
2(1 — nm) — (1 — 2n) : 
Of Go 
Gers 1 +(= — (1 <a 2n)T ex yi" —_ ra\(S)” a 
Ty Tr Go 
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A Method for Improving the Inherent 
Stability and Control Characteristics 
of Helicopters 


R. H. MILLER* 
Massachusetts Institute of Technology 


SUMMARY 


The problem of helicopter control and stability is examined 
with a view to establishing whether satisfactory inherent stability 
and control characteristics may be achieved without major design 
modifications and without having resources to automatic control 
devices. It is shown that the possibility does exist of improving 
both the damping and the static stability of a helicopter by a 
relatively minor modification of the blade mass and aerodynamic 
characteristics, together with the use of spring and viscous re- 
straint in the control system. This should result in considerably 
improved blind-flying characteristics and a reduction in excessive 
control sensitivity without sacrificing maneuverability. The 
control characteristic of such an inherently stabilized helicopter 
is evaluated by means of the transient response characteristics 


to abrupt control manipulation 
Section (I) 


INTRODUCTION 


5 es CONTROL AND STABILITY CHARACTERISTICS in 
hovering flight of several different types of heli- 
copters have been investigated in reference 1, where it 
was shown that the handling qualities could be consid- 
erably improved by incorporating in the control system 
a suitably designed automatic pilot. It was further 
suggested that satisfactory control could be achieved 
without the use of an autopilot by overbalancing the 
blade about the feathering axis and incorporating spring 
and viscous restraint in the control system. That the 
blade chordwise balance condition has a pronounced 
effect on the control characteristics has been known for 
some time from flight experience. Its effect on the 
damping in pitch of the helicopter was first demon- 
strated theoretically in reference 2. Its effect on the 
complete stability and control characteristics of the 
helicopter was first analyzed in reference 1, and the 
effects of viscous, aS well as spring, restraint were in- 
cluded. The present paper presents a continuation of 
this investigation and shows that what might be con- 
sidered ideal control characteristics may be obtained 
by this means, together with some further modifica- 
tions. The investigation is also extended to cover, in 
an approximate manner, the forward flight régime. 
Presented at the Rotating Wing Aircraft Session, Seventeenth 
Annual Meeting, I.A.S., New York, January 24-27, 1949. Re- 
vised and received January 5, 1950. ‘ 
* Associate Professor of Aeronautical Engineering 


The advantages of the proposed system of stabiliza- 
tion lie in its simplicity and also in the fact that a dif- 
ferent degree of stabilization may be achieved in pitch 
and roll, since it is proposed to apply the spring and 
viscous restraints to the nonrotating elements of the 
control system. It should thus be possible to stabilize 
the ship adequately in pitch without overstabilizing 


in roll. 


DISCUSSION 


Probably the most difficult part of any investigation 
of helicopter stabilization is the necessity of establish- 
ing what could be considered ideal control characteris- 
tics for which to aim. In the full realization that the 
definition is controversial, it is proposed to establish 
as a criterion the fundamental consideration that the 
helicopter motion should follow the control motion 
with as little distortion as possible. 

Taking for illustration, for reasons that will be dis- 
cussed later, the angular displacement in pitch of the 
hovering helicopter in response to an abrupt control 
displacement (unit step function) followed by neutral- 
ization of the controls 1 sec. later, this means that the 
angle of inclination of the helicopter must always be 
proportional to the displacement of the control stick 
from the neutral position, and this displacement must 
be achieved rapidly (for rapid control motion as as- 
sumed here) with a negligible overshoot. Furthermore, 
when the controls are neutralized rapidly, the helicopter 
must return rapidly to its initial trim position. Since 
the forward velocity of a helicopter is proportional to 
the angle of inclination of the tip path plane relative 
to the flight path, such a helicopter, when initially 
trimmed for hovering flight, would return to a steady 
hovering condition or any other trim position upon 
neutralization of the controls. Although this may 
appear to be a severe requirement in view of present 
helicopter handling qualities, it corresponds to the 
initial control characteristics of most conventional air- 
planes with respect to their flight path, as shown in 
reference 1 from which Fig. 1 has been taken. This 
may be compared with the hovering control character- 
istics of a conventional single-rotor unstabilized heli- 
copter as shown in Fig. 2, which has been taken from 


reference 3. The major causes of the present unsatis- 
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factory handling qualities are clearly indicated by these 
curves and arise from (a) the overshoot following 
neutralization of the controls, (b) the existence of an 
unstable oscillation whose period is sufficiently long 
to be manageable but sufficiently short to effect the 
initial control response, and (c) the existence of a large 
lateral motion following forward stick deflection. 


As pointed out in reference 1, the existence of an un- 
stable oscillation with a period of about 10 sec. is not 
of itself important, since few pilots would wait until 
the oscillatory mode established itself before applying 
corrective controls. However, the possibility of such 
an oscillation occurring is a clear indication of unsatis- 
factory initial response characteristics, since the condi- 
tions that permit the existence of such an oscillation are 
those that primarily effect the initial response and, 
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hence, the pilot’s opinion of the handling qualities 
This will be evident from an examination of the stryc. 
ture of the transient response equations given in Sec. 
tion (II). 

In reference 1, it was shown that these control char. 
acteristics may be considerably improved to the point 
where they could be considered almost satisfactory 
according to the definition given above if the damping 
in pitch of the helicopter is increased, and additional 
static stability is provided, by means of a suitably de. 
signed automatic pilot. However, the necessity of pro. 
viding such a device on all helicopters would result in a 
severe weight and cost penalty, particularly on small 
ships, and additional servicing problems. In attempt. 
ing to improve the inherent stability characteristics 
to the point where the automatic pilot could be elimi- 
nated, several different rotor configurations were con- 
sidered, including rigid rotors, rotors with a large 
flapping hinge offset, restraint about the flapping hinge, 
and the effects of modifying the aerodynamic and mass 
characteristics of the blade, together with restraint 
about the feathering hinge. This last method proved 
the most effective, and it was found possible to increase 
the damping in pitch by overbalancing the blade in a 
chordwise direction and providing elastic restraint 
about the feathering axis. The degree of damping 
realized with elastic restraint only was of the same order 
of magnitude as could be achieved by offsetting the 
flapping hinge appreciably from the center of rotation 
but was not accompanied by the additional restoring 
moment due to forward velocity, and also the additional 
control power, associated with such an offset. Conse- 
quently, the period of the slow oscillation was increased 
considerably to the point where its effect on the initial 
control response would be negligible, and the oscilla- 
tion itself became slightly convergent. On the other 
hand, the control sensitivity was reduced to about 
one-fifth of its previous value. The addition of viscous 
restraint in the control system was therefore proposed 
in order to delay the damping action and improve the 
initial response, and this is investigated in detail in the 
present paper. With the addition of viscous restraint 
(Fig. 3), the control characteristics are almost satisfac- 
tory except for a slight transient oscillation and except 
for the fact that, as a result of the restoring effect 
of the rearward inclination of the tip path plane due 
to the forward velocity, the helicopter would tend to 
return to its initial trim position even if the control 
deflection were maintained, requiring additional trim 
for forward flight. As a further consequence, when the 
controls are neutralized, a certain amount of overshoot 
would develop before the helicopter returned to its 
trim position. The next step is then to eliminate this 
restoring force with respect to forward speed by means 
of camber in the blades. This results in the control 
response of Fig. 4 which could be considered satisfac- 
tory except for the transient oscillation, which is still 


present. Physically, the stabilization arrived at by the 
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Fic. 5b (bottom). 
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methods discussed so far has been achieved by increas- 
ing the damping of the tip path plane relative to its 
ial plane of rotation, thus preventing it from follow- 
he control motion, or the helicopter motion, as 
This, of 


init 
ing t : 
rapidly as in a conventional helicopter. 
course, reduces the control sensitivity, as has been 
n, but provides damping moments on the helicopter 


see 
and, hence, 


itself, since, whenever the tip path plane 
the thrust—is not lined up with the helicopter, moments 
exist on the helicopter. However, when the damping 
of the rotor relative to its initial plane of rotation is 
increased, its damping relative to the helicopter is de- 
creased. Consequently, the rotor follows the helicopter 
in an oscillatory manner, giving rise to the high-fre- 
quency transient response shown here. 

The next step would then logically be to increase the 
damping of the rotor relative to the helicopter, and this 
may be done by introducing a feathering action on the 
rotor proportional to the rate of change of the angle 
between the tip path plane and the helicopter. In addi- 
tion, an elastic restraint may also be provided as a 
referencing device. Such a linkage physically would 
consist essentially of a link containing a spring and a 
damper (see Fig. 6) connecting either side of the swash 
plate (depending on whether positive or negative 
damping is desired) to a point on the blade outboard of 
the flapping hinge. The restraint about the blade 
flapping hinge provided thereby would be negligible, 
and its direct effect on the flapping of the blade may 
therefore be neglected. By a suitable choice of all the 
various stabilizing devices discussed so far, it is then 
possible to achieve a control response as indicated in 
Some slight distortion is still evident but is 
In addition, the 


Fig. 5. 
not believed to be of any importance. 
control sensitivity may now be increased to any desired 
value without distorting this response (Fig. 5b). 

No attempt has been made here to investigate the 
whole field, and it is possible that some other combina- 
tion of the various parameters discussed in this paper 
may be preferable. The powerful methods developed 
for servo synthesis have not as yet been applied to this 
problem, except in an elementary fashion in reference 1, 
mainly because the stability determinant for the heli- 
copter alone contains too many cross coupling terms 
(feedback) to permit a ready clarification of the problem 
by this method. 

Some justification for the method of analysis used 
here and in previous papers on this subject should be 
presented. As an illustration of the helicopter control 
response, the angle of inclination of the helicopter shaft 
axis—and, hence, of the pilot—has been chosen as the 
representative mode in establishing the transient re- 
sponse, since in hovering flight the pilot’s opinion as to 
the helicopter control characteristics would probably 
be closely associated with the motion of the horizon, or 
any fixed reference point, while attempting to hover in 
Also, a unit step function has been 
Any other mode, such 


a fixed position. 
assumed as the control input. 
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AXIS j TO PILOT'S CONTROL 


NOTE: 1) £ ~ (a, = 4.) WITH 90° LEAD OR LAG FOR + ke AND ky 


2) PILOTS CONTROL MUST BE RESTRAINED OR IRREVERSIBLE 
TO ENSURE STABILIZATION BOTH STICK FIXED AND FREE 











Fic. 6. Schematic sketch showing possible control linkage 


as the helicopter translation, may be used with any 
other type of input, such as an exponential control dis- 
placement (which would probably be more representa- 
tive of an actual pilot-applied control displacement), 
but it is believed that, pending further investigation, 
this is the most convenient method of synthesizing the 
complete system. Furthermore, once the response to a 
unit step function has been established, the response 
to any other type of input may be readily established 
graphically by the use of Duhamel's integral.‘ The 
velocity and acceleration of the helicopter may be 
readily obtained, if desired, by differentiation of the 
displacement response equation given in Section (II). 
This may be done more conveniently than the reverse 
operation of obtaining displacement by integrating an 
acceleration response. 

Although the problem of forward flight has been 
briefly considered in this paper, the investigation is 
limited chiefly to the case of hovering flight for the fol- 
lowing reasons: 

(a) Except for the problem of static stability with 
respect to forward speed as discussed in reference 5, the 
dynamic stability characteristics of the helicopter are 
not greatly affected, within the limits of present aero- 
dynamic rotor theory, by forward flight. The princi- 
pal change arises from the instability of the rotor with 
respect to its angle of inclination relative to the flight 
path, which in effect simply increases the instability 
of the whole helicopter to the point where the unstable 
oscillations may become simple divergences. A method 
of allowing for this effect is given in Section (II) of this 


paper. Further effects arise from the change in rotor 
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thrust with its angle of inclination to the flight path, 
which may be important when the helicopter center of 
gravity is displaced from the shaft axis. However, it 
would appear difficult to evaluate correctly such changes 
in thrust without introducing the rotor speed, or pitch 
setting if constant power and rotor speed are assumed, 
into the equations of motion. 

(b) The effects of the rotor slipstream on the fuse- 
lage at small advance ratios are not subject to rigorous 
analysis and have a pronounced effect on the stability 
characteristics. 

(c) The conventional theory for the induced flow 
through the rotor is completely inadequate for handling 
the nonstationary problems arising in forward flight 
stability analyses. 

(d) It is believed that the condition of hovering 
flight presents the most serious control problem, since 
it is a flight condition peculiar to the helicopter and one 
whose great potential advantages in blind flying cannot 
be fully realized until it may be used as a basis of orien- 
tation by the pilot. In forward flight the helicopter 
may be referenced to the flight path and the control 
position referenced to the forward velocity by means of 
tail surfaces or by means of the promising development 
reported in reference 6, which eliminates the problems 
of the slipstream immersing the tail surface at low for- 
ward speeds, or by means of blade camber as discussed 
in Section (II) of this report. In hovering flight, how- 
ever, the helicopter can only be referenced to earth 
axes if it is statically stable relative to earth axes, and 
this is the type of stability which it is believed must be 
provided. In practice, this means that, when in doubt, 
the pilot can release the controls with the assurance 
that the helicopter will return to a stationary hovering 
condition rapidly and with only such translation 
motion as may arise from wind velocities relative to the 
ground. 

The results of the analysis presented in this paper 
have been, with the exception of Fig. 2, for dual-rotat- 
ing helicopters. This permits considerable mathe- 
matical simplification with little loss in generality. 


CONCLUSIONS 


This paper has shown that the possibility exists, at 
least theoretically, of stabilizing a helicopter to any 
desired extent by relatively simple modifications to the 
control system and the blades. It is believed that a 
sufficient number of variables have been provided in this 
analysis to obtain any desired control response for any 
flight condition on conventional-type helicopters. 

However, the investigation reported here can only be 
considered as an introduction to an extremely complex 
problem. It is not presented as a final solution to the 
question of helicopter stability and control and, at best, 
is simply a theoretical indication of what might be a 
promising field of flight investigation. Modification of 
the blade balance characteristics will undoubtedly in- 
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troduce new vibration and flutter problems and accenty. 
ate the old ones, although many methods of minimizing 
these undesirable vibration characteristics may be de. 
vised and have been discussed in Section (II). Also, a 
wrong choice of the parameters discussed here and ¢e. 
rived in the theoretical portion of this paper could rp. 
sult in a decrease rather than an improvement in the 
inherent stability of the helicopter. One example of 
this has been given in reference 3 where it is shown that, 
when the flapping hinge eccentricity is large and the 
mass unbalance incorrectly distributed along the blade 
length, an unstable condition could occur. Further 
considerations arise from the incorrect use of blade 
camber, which may result in an unstable stick position 
with forward speed, from the incorrect use of an offset 
aerodynamic axis and camber that would result in ay 
unstable rotor with respect to its rotational speed and 
from the elastic problems of blade twist in the steady 
state condition and bending in the plane of rotation 
(discussed in reference 1). Also, suitable stops or a 
nonlinear spring restraint would have to be incorporated 
to prevent excessive feathering of the blade under 
flight conditions that may produce blade stall. 

For these reasons, in order to avoid incorrect gen- 
eralizations, no attempt is made here to give quantita- 
tive results for any specific problem. Methods are 
given in Section (II) for modifying conventional hel- 
copter blades to obtain the results discussed here and 
for calculating the desired degree of stabilization. 
These formulas should be used for a rigorous analysis, 
taking into account the peculiarities of any particular 
design, before an attempt is made to stabilize the heli- 
copter by the methods discussed in this paper. 


Section (I1)—Theoretical Development 


EQUATIONS OF MOTION 


In reference 1, the equations of motion for several 
different types of helicopters have been developed in 
considerable detail, the only assumptions used being 
(a) that the blades were infinitely stiff in bending and 
torsion, a necessary assumption for a precise solution; 
(b) that the inflow was constant over the rotor disc; 
and (c) that the aerodynamic derivatives could be ob- 
tained, as has been customary in airplane stability 
analyses, from stationary flow theory. 

It was then shown that the complete solution could 
be simplified without any appreciable loss in accuracy 
by neglecting derivatives arising from the accelerations 
of the tip path plane in pitch and roll and also the 
velocity of the tip path plane in roll when it is desired 
to examine the stability characteristics in pitch. Fur- 
ther simplifications, permissible when the flapping hinge 
is located close to the center of rotation of the rotor, as 
is usually the case, were also indicated. 

Since (tip path plane roll) is no longer a free vati- 
able, Eq. (17) of reference 1 could be written 
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a CBo ay 
HLH, — Hy yu - 
0 A oil 
Thay + TT(ay — B,) = 1750; (1) 
CBo ay 
(sm - Mo A )s + Q? + 
(Mt + Wy) (ay _ By) salon M541 (2) 
26 
Du ~— A(ay — B,) + = —Ab; (33) 


Q 


These equations of motion were set up in such a way 
as to allow for the fact that the thrust is not always 
perpendicular to the tip path plane, the primary effect 
of this phenomenon being to reduce the damping in 
pitch (or roll) of the helicopter. (See reference 1, page 
465.) A further effect is to add horizontal forces in the 
tip path plane due to blade coning, Bo, although there 
may be some question as to the precise magnitude of this 
increase when the blade is extremely flexible in bending. 
This latter effect, however, is not of prime importance 
when considering the dynamic stability and the control 
characteristics, although it is of importance when 
considering the static stability problem of stick position 
versus speed. The reduction in damping, however, 
becomes important when there is no offset of the flap- 
ping hinge (teetering or rocking rotor) and when the 
rotor blades are designed to operate at a relatively low 
angle of attack (high solidity, high tip speed, light disc 
loading). 

The left-hand side of Eqs. (1) to (3) may be expanded 
into the following stability equation when the harmonic 
substitution e” is used: 


Ay = Cyt + Cy® + Cv? + Cw + Cs (4) 
where, to a sufficient degree of accuracy, 
C, = 2, C, =A 
C; = (MM + 9G)(2 + nD) + 


CBo C D 
A | »( 20: — Wy * +h~ + ; 1 | 


a a 


C, = 0 
Cs = Hof (M + M)D + Al[IM, — Mo(CB/A)]! 


The nomenclature is that of reference 1. 

With the same assumptions as were used in deriving 
Egs. (1) to (3), the equation of motion for an infinitely 
rigid flapping blade about the feathering axis, when the 
blade chordwise c.g. and the aerodynamic center are 
not coincident with the feathering axis and when the 
blade pitching moment coefficient, Cm,, is not zero, may 
be written* 


Ky Q 26; 
7,02”! + [,Q? 6, = o —A (a, — Bi) + 
b2R3 ‘i 
u(p + os Cno) + Ad; (5) 
2Tr 
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where A, and ¢ are the elastic and viscous restraints in 
the control system. 

The inertia of the blade about the feathering axis has 
been neglected, since it was shown in reference 4 that 
this approximation is justified, at least for the case 
where the natural frequency of the nonrotating blade 
about the feathering hinge is greater than the rotational 
speed of the rotor, which would usually be the case. 

Substituting Eq. (3) in Eq. (5) and writing 

ky = (27,2/K,)[1 — (A’/A)] 
Ky A 2Tp 
h = a/Ky 


Il 


results in the control equation 
+ hb = kB, + Ro (6) 


It is thus evident that the effect of offsetting the aero- 
dynamic and gravity axes from the feathering axis, 
providing camber, and restraining the blade about the 
feathering axis is equivalent to introducing an auto- 
matic pilot in the control system with a characteristic 
time lag, /;, and with an output proportional to the 
rate of pitching of the tip path plane, 8,, and propor- 
tional to the rate of displacement of the helicopter, x, 
from its initial hovering position. 

The predominant term in 2 is that due to Cm, since 
the other terms arise from the usually negligible, first 
harmonic sine component of the blade lift. 

If the viscous restraint, ©, is provided in the non- 
rotating portion of the control system (for example, from 
the stationary rim of the swash plate to the helicopter 
frame), then Eq. (6) may be written in the algebraic 
form 


(Oa) 


6, = (RvB, + Roy) (1 + 1X) 


where \ = vQ in the harmonic substitution e””. 

The above form of Eq. (6) is for the case of springs 
and dampers in parallel. From the servo point of view, 
more satisfactory stabilization could be obtained by 
means of springs and dampers in series. Assuming ke 
to be zero, this arrangement would give an equation of 
the form 


A= kB; + kp; 


which is similar to the automatic pilot considered in 
reference 1. This system has not been investigated 
thoroughly because it is not considered desirable to 
transmit. control forces directly through a viscous re- 
straint without incorporating an elastic restraint in 
parallel because of the possibilities of failure and of 
severe nonlinearity of the viscous element and also 
because, if k2 were not zero, then the possibility of a 
slight divergence exists when ke is made negative. 

If the swash plate is connected, through a spring A» 
and damper ¢ in parallel, to the blade at such a point 
as to respond to the blade flapping displacement rela- 
tive to the helicopter (a; — §;), but in such a manner 
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that the pilot's control does not work against these where A, is given by Eq. (4) and, to a suflicient degre 
springs (such an arrangement is indicated schematically of accuracy 
in Fig. 6), then Eq. (6) may be written 


a = Do. 
A + 16, = RsBi + Rear — Bi) + b = A(M + Mo) 
Ry( au 7 Bi) + Regu (7) CBo 
4 . = ° . Cc = A(t + me) (HH aa Hy, - ) — 
where ks, ke, kz, ks, and /, are functions of Ky, Ko, ¢1, C2, A 


and the geometry of the control linkages. Their 
derivation for any chosen arrangement presents no 


CBo 
AH, (on. =e Jt | —_ DH sm 


particular problem. It is, however, important that the p =A 

viscous restraints, ¢, and ¢, be linear. The friction in q = 0 (may be assumed) 

conventional dashpot-type dampers obviously detracts r = A(M + Mo) 

from their linearity, and resource would have to be 5 = A(M + Mo) [Ly — Hy(CBo/A) — (D/A) 
made to devices similar to those developed at the w, y = 0 (may be assumed) 

M.1.T. Instruments Laboratory, where such integrating x = AM) + An(M + M) 


means have been developed to a degree of perfection 
far beyond that required for this application. 

From Eqs. (1), (2), and (3), it is possible to derive 
the following relationships: 


= AH,(M + No) 


a 


If the pilot moves the control stick forward and thereby 
applies a cyclical control deflection, @,, to the helicopter 
control system through a linkage system similar to that 


a/—0, = (av? + by + c)/Ay (5) of Fig. 6, then, from Eq. (7), 
Bi/—6, = (pv? + qv? + rv + 8)/Ai (9) NksBr he + kz ks 

= @, — (a, — pi) + (7 
u/—O, = (wv? + xv? + yo + 2)/ A) (10) i+—”' 14+K™ "1+ 9 


Eqs. (8), (9), and (10) may be substituted in Eq. (7a) to obtain the result 


my (1 + Ld)(av? + bv + ©) , 
—0, Al + led) + [(hs — By)A — Bel (pv? + rv + s) + (Re + Abr) (av? + bv +c) + Ae(xv? +2) * 


Eliminating v from this expression by means of the relationship A = vQ (in order to permit the control respons 
to be expressed in seconds) results in 
a,/—0. = [(1 + brd)(ad?2 + BOA + cN?)O4 / Ay 
where 
As = BB,» + Bov4 + Bs? + By? + BsX + Bg (Lla 
and, to a sufficient degree of accuracy for most applications, 


B, = Cie 


By = Cy + 1QC2 + (ks — bz) Qp 
Bz = CQ + 122°C; — RQp + kQ*a 
By = CQ? + RgrQ® + ReaQ? + xh, 
Bs = 1.2'C; + (ks — z)sQ* + RycO4 
Bg => C;24 + ReQ4(c — §) + R24 


In establishing the initial transient response of a helicopter (over a period of time of the order of 10 sec.), tht 
terms a and c in the numerator of the right-hand side of Eq. (11) may be neglected. 

The fifth-order equation, represented by As, has five roots, of which one root represents a convergence that 1s 
usually sufficiently rapid to be negligible in establishing the transient response. When this is the case, the larg 
root is, approximately, 


As = — Bz, B, 


and need not be considered when calculating the transient response. Under these conditions Eq. (11) takes tle 


form 


~~ - B; Be (NEMA + a)(A F Ag(A + Ad) (12 
4 3 2 . 
eS we Sey 5 +e 


Qa (=) a Isr + l : lor + ] 
—0,\b93) 
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It should be emphasized that this form of Eq. 11 is a rough approximation and should be used only when \s is 


rea: 


for t 
simplification of the transient solution. 


sonably approximated by the first two terms of A,. Consequently, \; to 4 should always be determined from 
the complete equation represented by As, by first extracting As, in order to check the validity of the approximation 
he particular case being considered. The use of this approximation, when justified, results in an appreciable 


TRANSIENT SOLUTIONS 


The roots \; to \4 are usually complex and may be written 


Mi == A + 1B, 


de = A —_ 7B, 


As = P + iQ, \4 = = iQ 


The transient response, when @, is in the form of a unit step function, may then be obtained, by the usual methods 


of operational calculus, in the form 


ay (on) Wh, as (F ;Ot + # xs m) - (F ; Bt A. Bt) 
—9.\0Q3 wae” cos ¢ 0 sin ¢ G cos B sin 


where 


F = 2(P — A) —1(P? — 


G = (P — A)'+ 2(0? + B*)(P — A)? + (@ -- B?)! 


A?) — 1,(Q? — B?) 


J =(P — A)? — (@? — B’) — 24Q0°(P — A) — hP(P — A)? + 12P(Q? — B?) 
K = (P — A)? + (? — B’) + 24.B°(P — A) — 1A(P — A)? — hA(Q? — B?) 


Some duplication in nomenclature occurs as a result of the above choice of symbols, but it should not cause any 


confusion. 


Several transient responses have been obtained by use of this equation, one of which is shown in Fig. 3. 


When the long-period oscillation is eliminated by use of sufficient negative ks to make \, = 0, then the three 


remaining roots may be written in the form 
Ay => A + iB, de 
and the transient response becomes 


a (=) (4C — 1)e-* 
4 hee "th 
-6,\003/ C(A® + B?) 


e Al 


[B(C = 2A)]|? + [AC — (47 — B?)}? 7 


C[B?+(C— A) [B?+ 


A-—-1B, 3=C 
le-* sa. 
cos Bt — sin Bt} — 
(C — A)? B / 
AC — (A? — B? 
Yic — 2A) cos Bt + | B J sin Bit 


Two examples of transient responses which may be obtained with this equation have been given in Figs. 4 and 5. 


These, together with Fig. 3, have been discussed in Section (I) of this paper. 


FORWARD FLIGHT 


In steady level forward flight, the main differences 
in the equations of motion arise from the fact that the 
rotor becomes unstable with reference to angle of 
attack, as pointed out in reference 5. A rigorous analy- 
sis of the forward flight régime should include: 

(a) The degree of freedom represented by the rotor 
speed with constant engine power or, if a constant speed 
rotor were feasible, the change in collective pitch, 4, 
with changes in rotor inflow. 

(b) The resultant changes in the induced flow 
through the rotor. 

(c) The vertical degree of freedom of the helicopter. 

(d) The degree of freedom represented by the coning 
angle of the rotor, which is no longer constant. 

(e) The influence of the fuselage, including slip- 
stream effects. 

Such an analysis would be laborious and _ scarcely 


justified, in view of the present limitations in rotor 
aerodynamic theory, unless wind-tunnel data on the 
stability derivatives were available. However, as a 
rough approximation and since it is desired to discuss 
only the initial transient response, the steady-state 
conditions may be assumed to persist, and the only 
effect that need be considered is the rotor instability 
with respect to angle of attack. This may be approxi- 
mated by adding a term —2A oA; to the left-hand side 
of Eq. (3), which in turn changes C, in Eq. (4) from 
zero to 

Cy = —2A(M + My) uo? 
and adds a term —2A(IW + 9My)uo?Q2* to B;, and the 
same term, multiplied by /2, to By in Eq. (lla). yo is 
the steady-state advance ratio. Evidently additional 
k; is all that is required to reduce the equations to their 


original numerical value for any given case. Additional 


kg would also be needed, when it is desired to eliminate 
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the long-period oscillation, because of the effects of 
rotor and fuselage drag when yo is not zero. 


Eq. (5) is also unchanged except for the addition of a 
term 


yo(pb2R®C,,/2Tp) 


to the right-hand side. 


Evidently, negative C,,, will result in an appreciable 
reduction of stick motion with forward speed which may 
reduce to a condition of instability with forward 
speed. In general, the choice of the stabilizing 
parameters would have to be the result of compromise 
between the hovering and cruising stability characteris- 
tics. 


DESIGN PROBLEMS 


When mass or aerodynamic unbalance is applied to a 
blade, the vibration problems of the helicopter may be 
expected to become more serious. In particular, the 
third harmonic oscillations, which are critical on a three- 
bladed rotor and which arise primarily from second 
harmonic disturbances, would be more severe, since a 
second harmonic feathering action of the blade would 
occur if the blades were flexible in torsion. However, 
if the blades were stiff in torsion, the dampers c, and 
c, would effectively cancel these higher harmonic oscil- 
lations. These dampers will also minimize the second 
harmonic oscillations in the fixed system (swash plate) 
because of first harmonic feathering of the blade in the 
rotating system, which would be critical on a two- 
bladed rotor. 


Blade flexibility in bending will also considerably 
modify the term k; appearing in Eq. (6). In particular, 
for the hypothetical case of an infinitely flexible blade 
with no curvature (mass and aerodynamic forces equal 
and opposite at every point along the blade), the term 
A’/A would be unity (and k, then zero) when x, = x;. 
This is not the case for a conventional blade with ap- 
preciable stiffness or even for a highly flexible blade 
since an elastic curvature of the blade would always 
exist. 


* These considerations of blade flexibility and others, 
such as the necessity of canceling out the steady twist- 
ing moments and the effects of blade curvature in the 
plane of rotation, as discussed in reference 1, would 
indicate the advisability, in the case of large flexible 
rotors, of using a small “‘pilot rotor’ for control and 
stabilization such as was proposed by Bruel for the 
French helicopter SNAC Nord 1700. 


In order to permit a wide variation in the parameters 
of Eq. (6) or (7) and also to permit ready modifications 
as indicated by flight experience, it would appear logical 
to obtain the parameters &,, ke, and /; by means of sepa- 
rate masses and auxiliary aerodynamic surfaces located 
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at variable spanwise distances from the flapping hinge 
and variable chordwise distances from the feathering 
axis. 

If a mass 1, is located at a radial distance R, from the 
flapping hinge and a chordwise distance x, from the 
feathering axis, then 


Ip = m,x,R, 


If an auxiliary surface of span d, and chord by, wit, 
a lift curve slope a4, is located at a radial distance R, 
from the flapping hinge, has a chordwise distance froy 
its aerodynamic center to the feathering axis of Xs 
and is set at an incidence 4, relative to the main blade 


Pe pa4b,R* (=3*) (5) 
Tr R' 2 


— pa,b,R* (ax Ry Ar {I 
D —_ a ) R? | + 604) R —_ 9 { 


Cmo4 = 14004X4 bs 


then 


and the expression containing C,,, in Eqs. (5) and (jj 
becomes 


1 /bR*\ pa4b,R* (daxaR, 
pe p c a 3 Bo. 
“ Tr Tr R 


where FR is the rotor radius. 
when ahead of the feathering axis. 


xX, and x4 are positive 


In choosing the variables in the above expression in 
order to obtain the desired stability and response char 
acteristics, it should be realized that: 

(1) Too much negative 4), will result in instabilit 
of stick position with forward speed. 

(2) Negative A’ in conjunction with positive 6, 0 
positive A’ with negative 0), may increase the difficulty 
in maintaining constant Q. 

(3) The torsional stiffness of the blades should bx 
high. Alternatively, the steady-state moments mus 
be canceled out on the blade itself,* which in tum 
would thrust with for 
ward speed, since the coning angle 8» is a functio 


modify the variation of 
of po. 

In order to achieve both stick-fixed and stick-fre 
stability, the springs and dampers K, and ¢ should 
located (Fig. 6) in series with the controls and betwee! 
the swash plate and the helicopter frame. While tht 
dampers c; would undoubtedly eliminate much of th 
higher harmonic stick motions resulting from xa, x 


and @,, it may be desirable on large rotors to conttl f 


the rotor indirectly by means of ailerons or servo tabi 
located on the blades, as in the case of the Landgral 
and Kaman helicopters. 

Finally, it should be noted that small rotors with 
relatively stiff blades may be more easily stabilized }) 
the methods discussed here than the large rotors wit! 
blades flexible in torsion and bending. 
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CONTROL 


YOUNG STABILIZER BAR 


The Young stabilizer bar has been described and 
analyzed in reference 7. In order to compare its stabil- 
izing action with the system discussed in this paper, the 
following brief analysis has been made, with the simpli- 
fying assumption that the application is to a dual rotor 
helicopter. 

If By is the angle, at any azimuth position y, of the 
stabilizer bar relative to a system of axes fixed in space, 
if a, is the angle of inclination of the helicopter shaft 
axis at any position y, and if c, is the effective restraint 
between the helicopter frame and the stabilizer bar 
provided by the viscous dampers, then the equation of 
equilibrium for the bar about its rocking hinge becomes 


M.— M,=0 


i 


where ./, is the moment due to the dampers and M, is 
the inertia moment. 

If /, is the moment of inertia of the stabilizer bar 
about its rocking axis, then 


M, = T,Bsy, M, = Clay — By) 


If, as was assumed in deriving the equations of motion 
for the blades, 


Bsy i Bs, cos vy +- Bs sin y, 


a, = a cos yp 
and writing 


as, = Qj — Bs b= Byoy c= 6, i. 


SI 


the equations of motion become, after reducing to alge- 
braic form by means of the harmonic substitution e”’, 


s,(Qe + 2Or) + bs,(Ac + A) = 2Odaxy 
as,(Ac + A®) — b,,(Qe + 2Od) = AZay 


whence 


ds, _ 20°A(c + 2X) + AKC + A) 


a  2%(c + 2A)? + AX(c + A)? 

It is generally of interest to determine the effect of 
the stabilizer bar on the small roots only of the stability 
equation for the complete helicopter. Consequently, 
would be small compared to 2, and the above expression 
could then be approximated by 


Qs, 2r = (2/c)Xr 
a ctHt2r 14+ (2/c)r 


Since the relative motion, as;, between the helicopter 
frame and the stabilizer bar is used to apply a cyclical 
pitch variation to the rotor and since, by means of an 
ingenious linkage arrangement, the pilot-controlled 
feathering action, 6,, is made independent of this cycli- 
cal pitch variation, it follows that Eq. (6a) or (7a) 
would become, for this case, 


A, = 0. + [Rida /(1 + AA)| 


where k; is a simple function of /; and depends on the 


link ratios and c. It is thus apparent that the Young bar 
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produces stabilizing characteristics similar to those dis- 
cussed in this paper in connection with Eq. (6), except 
that ky = 0, ki = f(i,) and is proportional to a and 
not Bh. 


HILLER CONTROL ROTOR 


The Hiller control rotor has been described and 
analyzed in reference 8. 

Using subscript 7 to denote the control rotor param- 
eters, Eq. (3) for the control rotor, when applied to a 
dual-rotating rotor, becomes 

9 
— Aa, — Bi) + oY os —A,6, 


Q 


Diu (3a) 


Using the definitions of X and P of reference (2) and if 
6, is the feathering action applied to the rotor by the 
pilot 


6, = 0. + X(Bi — ay) 
and 6; in Eq. (3) becomes 
6, = P(a, — B,,) 
whence 
Bi, = & — (6; P) 
and Eq. (3a) may be written, in terms of 6;, as 
2 2P D 


AQ ian r A, 


r 


— 

PX (a —_ 2) + P6. 
Writing 

ls = 2 AQ, 

kg = P(D./A,), 


ks == 2P AQ 
ke = —PX 


This equation may be expressed in a form equivalent 
to that of Eq. (7a) as 


on ( P ) ee 
eee "Bee ee 
k 


7 ks 
1 + Jd 


— 1) 
(ay UF pat 


and it is apparent that the control rotor performs in a 
manner somewhat similar to the stabilizing devices 
considered here, except that k; = 0 and &; is propor- 
tional, as in the case of the Young stabilizer bar, to 
a, and not Bi. However, the pilot applied control, @,, 
is subject to the time lag, /:; this is not the case for 
the Young bar and the stabilizing devices considered 


here. 


KAMAN SERVO CONTROL 


The Kaman servo control has been described in recent 
articles, but no analysis of the apparent stabilizing 
effects of this control system is, to the writer’s knowl- 
edge, available. However, because of the location of 
the servo tab aft of the blade elastic axis and in view of 
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the elastic restraint provided by the blade torsional 
flexibility (which is the only feathering freedom pro- 
vided), it is probable that a certain amount of the 
parameter k, in Eq. (6) with relatively little /; is inher- 
ent in the rotor system. It is therefore possible that the 
Kaman helicopter has a greater amount of damping in 
pitch than most conventional helicopters, with the 
consequent improvements in the stability and control 
characteristics. 

Since the Kaman helicopters employ counterrotating 
rotors, the analyses contained in this paper would 
probably be more applicable to this design than to the 
Bell and Hiller machines. In the latter two cases, it 
is to be expected that the stability characteristics 
would be appreciably modified by the coupling between 
the pitching and rolling motions of the rotor tip path 
plane. 


NUMERICAL COMPUTATIONS 


Figs. 2 to 5 have been calculated for a helicopter with 
the same parameters as were assumed in reference 1. 
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(For a list of these parameters see reference 
463 and 466.) 
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RIEF REPORTS of investigations in the aeronautical sciences 

and discussions of papers published in the JOURNAL are 

presented in this spec ial department. The publication will be com 

pleted approximately 8 weeks after receipt of the material. The Edi 

torial Committee does not hold itself responsible for the opintons 

expressed by the correspondents. Contributions should not exceed 
800 words in length 


Flow of Nonviscous Elastic Fluids 


F. G. Gravalos 

Department of Aeronautical Engineering, Rensselaer Polytechnic 
Institute, Troy, N.Y. 

March 27, 1950 


N THE ARTICLE “‘Steady Flow of Nonviscous Elastic Fluids in 

Axially Symmetric Channels”’ (JOURNAL OF THE AERONAUTICAL 
Sciences, Vol. 17, No. 3, p. 144, March, 1950), to obtain the 
equations of motion, M. H. Vavra uses the method of general 
orthogonal coordinates that was originally introduced as the 
natural system to study turboflow problems by the present writer 
as early as February, 1948 (see abstract, Bulletin of American 
Mathenatical Society, Vol. 54, No. 5, May, 1948). Mr. Vavra 
follows, up to his Eq. (26), conceptually the same procedure as 
that given in an O. N. R. Technical Report 3, Contract N7-ONR 
336-TO1, and distributed by the Office of Naval Research in 
August, 1948. 
and the rigorous application to the calculation of the flow through 
a turbomachine were presented in a colloquium given by the 
writer at the Massachusetts Institute of Technology in May, 
1948, a colloquium at which Mr. Vavra was present. 


The advantages of using orthogonal coordinates 


That the normal direction to the stream surfaces is a privileged 
direction, along which the flow may be integrated once the meri- 
dional flow is known, is the key analytical concept on which the 
solution of the problem is based; this concept is actually used 
whenever orthogonal coordinates are employed. This was also pre 
sented at the VIIth Congress of Applied Mechanics in London, 
September, 1948; later it was developed into a consistent theory 
published in Bulletin No. 62 of the Rensselaer Polytechnic In 
stitute, “A Laminar Theory of the Flow through a Turboma- 
chine.”". The technical publications mentioned above will be re- 
ferred to by T.R. (Technical Report), L.T. (Laminar Theory), 
and C.A. (abstract in the Transuctions of the Congress, Section 
II, page 599). 

Mr. Vavra obtains, by elementary methods—without using the 
concept of normal curvature—his basic Eq. (26) in the JouRNAL 
paper. This equation is a particular case of Eq. (47) of T.R., 
Eq. (59) of L.T., or of the equation in C.A. and may be obtained 
by setting w = 0 and e = 0 in the writer’s equations. However, 
the coefficients sin \/y and 1/R in Eq. (26) are, with the proper 
signs, normal curvatures of the stream surfaces. 

Neglecting the true meaning of the relation Vir = constant, Mr. 
Vavra proceeds to use his Eq. (30) to define the variation of V2 
along the normal to the stream surfaces by supposing that the 
Stream surfaces are known. 

In contradiction with Mr. Vavra’s interpretation of the prob- 
lem, for the case of an empty channel: 

(1) Vir = constant is an independent integral of the equa- 
tions of motion, U = 0 in Eq. (37) of L.T. 


(II) Eq. (26) is superfluous as an equation for V;, and his Eq 
(30), obtained from Eq. (26), is redundant as an equation for 
Vo. 

(III) For incompressible fluids, rotational or irrotational 
flows, the equation of continuity alone defines the meridional flow 
pattern (i.e., the meridional streamlines plus the velocity dis- 
tribution over them), and it does not contain V; implicitly or ex- 
Hence, V; should be determined from Vir = constant 
(Plus boundaries conditions. ) 


plicitly. 
and p from Euler’s integral. 

(IV) For rotational meridional flow patterns, whether the 
fluid is compressible or incompressible, none of the approximate 
methods quoted by Mr. Vavra possesses any applicability what- 
soever to the determination of the initial shape of the stream 
lines. 

Finally, let it be noticed that any theory of the flow through a 
wheel includes, as a special case, the flow through guiding vanes; if 
the wheel is supposed with infinitely many blades, then the corres pond- 
ing particular case is that of an empty channel. 


Approximate Incompressible Flow Solutions 
for Slender Polygons 


E. V. Laitone 

Associate Professor of Mechanical Engineering, University of 
California at Berkeley 

Received March 29,1950 


hes TRANSFORMATION FROM THE PHYSICAL 5 = xX + /y com- 
plex plane to the complex potential F = ® + 7i¥ plane for 
a uniform flow approaching a polygon symmetrical about the x 
axis, as shown in Fig. 1, is given by 

\ 
I 1 
a (F — ,)&/*) (1) 


n= 1 


| dz 
w dF U 


Since ¥ = 0 on the boundary and the flow velocity 1s ’ at x = 
+o, the transformation constant in Eq. (1) is 1/U as long as 
N 
z. (rf — o,) = VU 
n= 1 
that is, as long as the polygon profile is closed or concludes with a 
surface parallel to the free-stream direction as shown in Fig. 1. 
Now w is the conjugate of the complex velocity vector w = u + 
iv, SO 
dF 


\ 
w=u-w= =U ¢ 
dz 5 


(F = jl (on/r) (2) 


n 

The solution given by Eq. (2) is not directly applicable, since 
the relation between z and F must first be ascertained by inte- 
grating Eq. (1), which in most cases is impossible in terms of a 
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finite number of elementary functions. However, for a slender 
polygon wherein each (7 — a,) << 1, we may use the approxi- 


mate relations for the particular streamline ¥ = 0 
w=u-wru . 
= ; ; re (3) 
F= Ux+e+t+ iv = Ux{ 
so that Eq. (2) is simplified to 
N 
u oe 
_ e ey, | —ten/0) 5, nal (4) 
l ps l 
n= 1 


This provides a general solution at ¥ = ( for a slender polygon 
profile as long as 
N 


>» (rf — o,) = 0; 


n= 1 


5 — Gan\ << 1 
The results from Eq. (4) agree extremely closely with the 
available exact solutions for wedge or double wedge profiles, 
even when the value of (r — a) is as large as 7/4. They also are, 
for all practical purposes, identical to those given by the first- 
order solution for the two-dimensional source distribution repre- 

senting any profile by (see Sauer’) 

u 1 yo (EdE iw . 
~ l _— ~ ? (9) 


l 1 é—x l 


where the Cauchy principal value of the integral is to be taken 
at — = x and y’ is the local slope of the surface of the profile. 
For example, for the simple wedge shown in Fig. 2, we have, 


from Eq. (4), 


Ww 1 6/r 
as g8/™ |e — 1 |—8/* om 
l (L/x) - 1 
and, from Eq. (5), 
w ' $, L-x v4 adi 
cot n = n 
l 7 x (L/x) — 1 


These two results, identical in Fig. 2, are both within 1 per cent 
of the exact solution for any wedge having 6 < 1/10. Similar 
results within the same order of accuracy are also found for the 
double-wedge profile. 


REFERENCE 


' Sauer, Robert, Theoretische Einfiihrung in die Gasdynamik, p. 27, Eq 
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Remarks on the Velocity of Sound 


Max M. Munk 
U.S. Naval Ordnance Laboratory, Silver Spring, Md. 
March 22, 1950 


STRICT DEFINITION of the velocity of sound seems still to be 
in doubt. It seems to me that the velocity of sound c is 
defined by the necessity to give a Mach angle u leading to physi- 
I will now show that it does so 
provided the velocity of sound is taken in accordance with 


(A) 


the said partial differential quotient to be taken in the only 
uniquely defined direction 


cally significant Mach lines. 


c? = Op/dp 


namely, in the direction of the stream 
line of a steady flow. 

The proof requires that shear stresses be absent in the fluid. 
No potential or barotropic expansion law is required. The 
following proof is restricted to two-dimensional flows. The sym- 
bols have their usual meaning. 
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Fic. 2 (Laitone) 


At the point considered, let the velocity be in the Y directig: 
and let 6 denote the flow angle. 
the form 


The dynamic equations assuy 


(B) 
(C) 


uus + (p/p) = 0 
uv, + (py 'p) = () 
and the equation of continuity becomes 


(D) Upr + plz + ply = 0 





Multiply (C) by u, eliminate wu, from the system of equations 
(A) to (D); write the remaining equations 
(E) sin uw cos uw) = 0 
(F) (py sin w/c*p) + (0. cos pu 
and form (E) + (F) and (E) — (F). 
The equations thus obtained express the fact that along the 
Mach lines defined by (A) the characteristic differential 


(px cos u/c*p) + (Oy sin p 


sin uw cos uw) = 0 


(G) dp sin uw cos wh + dpc? = 0) 


is zero. No definition of the velocity of sound other than (A 


would lead to an equivalent relation. 


The Problem of Edges in the Doublet 
Distribution Method of Obtaining 
Supersonic Lift 


Theodore R. Goodman 
Cornell Aeronautical Laboratory, Inc., Buffalo, NY. 
February, 1950 


ECENTLY, Miss COHEN pointed out in reference | that the 
method of solving for lift solutions as given in references? 
and 3 “appears to be in error in the neighborhood of the leading 
edge or of a raked-out tip.”” The author is in agreement with 
this statement and intends to clear up the matter in this note 
In the appendix of reference 2 the solution to the integral equa 
tion 
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f(x) = 
is given as 
u(s) = 


This solution is not unique, however. 
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x—~-t 
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m the properties of linear equations it is clear that any solu 


Fro ; 
, will also satisfy Eq. (1). 


1 of Eq. (3), when added to Eq. (2) 


tior 


The solution of Eq. (3) 1s 
u(s) = A (s — a)! m (4 


which may be verified by substitution. Here, 4 is a constant. 


AN APPLICATION 


Consider the raked-out wing shown in Fig. 1, and apply the 
method of reference 2. Using the notation of reference 2, the 


integral equation to be solved is 
“! "uF “kv . 
1 di j f(my, v%)duy ACpi (11, v,)duy 
i F os 8/, 
| ele { } um My . (uu — ™) /? { 
. . 0 


Considering the 


0 
(oO 


Now apply Eqs. (2) and (4) to this equation. 
quantity in brackets to be the unknown, there is obtained 


a a "kn . , 
* f(a, m)du ACpi (1, 21) du A(u 
5 les 3], 

, == , _ ) oa 
J ku ul ul » (i uy Vv v 


Applying the solution again, there is obtained 


] ~ du 
fiz,4.) = — 
2x J; s— nu)’ 
JS ki ° 


2r VM kn \* 


"kv . 
AC (1, 0) du, 
3/s 
(u — uu) /* 
A(u)du B(v) 


(2— RV) 


Since the pressure coefficient cannot become infinite along the 











edge x, = 0, 4 = 0. Thus, Eq. (10) becomes, upon integrating, 
re hes “kv ACpi (1, 11) duy 
f(s,%) = — rs 
7 = 1 (f — mu) Vko — 1 
B(v) 
i (8S) 
2 — kv) 
tid -U, 
- yy, 
@,” 
7 
? 
v,=0 } 
@ ,” 
4 
4 
7 
ae v,zV 
£. —U,7U vv, 
hu zkv, 
Ulu 
Fog. 4. v 








The last term in this expression becomes infinite along the edge 
2— kv, = (0. In fact, for any region of any wing, the homogene- 
ous solution yields an arbitrary function of one of the Mach 
coordinates divided by the equation of the edge set equal to zero 
It is therefore necessary to examine 
the kinds of edges a region can have in order to see when the 


and raised to the '/. power 


arbitrary function is zero. 

There are three types of edges toa region: (1) a leading edge, 
(2) a trailing edge, and (3) an edge that borders another region. 
In type (2), if the Kutta condition is to apply, the arbitrary func- 
tion must be zero. In type (3), since the pressure coefficient can- 
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not become infinite on the interior of a wing, the arbitrary func- 
tion is again zero. For type (1), however, the pressure coefficient 
may become infinite, and it is necessary to have a method for 
evaluaiing the arbitrary function. 

It should be pointed out that the method of references 2 and 3 
is actually a method of acceleration potential. The boundary 
condition for 0a/Ox is satisfied rather than that for a; the Kutta 
condition is easily applied; and also there is an arbitrary function 
that must be evaluated from the condition on a, all properties of 
an acceleration potential. 

The arbitrary function arises because 0a /Ox is used [reference 3, 


], instead of a, as the basis for the analysis. The equation 


Eq. (7 
for a in terms of AC, can be derived from the following equation 
for the velocity potential given by Ward.‘ In the notation of 


reference 2, 


l 
g(x, y,2) = — x 
tr 
a 3 ACp(x1, Mi)s(x — x )dxidy, 
; (9 
Jr iy — mm)? + 27] V(x — 1)? — (y — yi)? — 2? 
Eq. (9) can be written 
1 Oo ; 
P(X, ¥, 5) = AC, (%1, 1) 
tr Oy Jr 
sv (x — x1)? — (y — m1)? — 2? 
tan”! dxidy,; (10 
(y — y1)(x —x) 
and 
Oo B oO 
OZ tr OY . ” 
(y — 1) (x —x1)[(x — 41)? — (y — 1)? — 227 AC, (m1, ni dxidy, 
(11) 
[(y — nn)? + 22) [(x — m1)? — 27IV (x — m1)? — (9 — nn)? — 3 
Now let Z ~ 0. Then 
() 
r = a= 
Z/ 7 ~ 9 
8B oO (x — x)* — (y — ni)? ; 
- \ = AC F (Xi, M \dxidy, ( 12) 
4m Oy r (x — m(y — n) 
By reasoning similar to that used in reference 2, one can 


derive, from Eq. (12), the integral equation 


re) . 
Af(xi, y)dxidy, = 
oy 2 
Oo - 
/ KAC, (x1, w)dxidy, (13) 
oy Ga 


V (x — x)? — (y— yn)? 


where 


(F — BHI — Hi) 


This is paralleled by the equation from which Eq. (5) was 
derived, except that it is based on an equation for a@ instead of 
Oa/Ox. 

If Eq. (8) is substituted into Eq. (13) [4C,: = constant 
4a/8], one arrives at an integral equation for B(™), 
B(v,)dry 


*y/(1 — k) 


=—— 3/e l1+k 
3/. y =a (14) 
411 — k)’4) Jo (2. -») Vi-k 
The solution is 
v / C 
B(z) = (l1—k) V1 +k Ve 4 (15) 
pr / 2 


where C is a constant. 
But C = 0, since AC, is finite along the line y, = 0. The com 


plete pressure coefficient in the tip region is, thus, 
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AC ACn — rr ky — “| 4 
2 = cos 2 
P wT V1 + uy 
ACh / v1 
(l1-k)V1+k¢ (16) 
a bd NV: — ky 


This is in agreement with the solution obtained by Evvard.® 


SUMMARY 


Three types of edges can border a region: (1) a leading edge, 
(2) a trailing edge, and (3) an edge bordering another region. 
For types (2) and (3) the method given in references 2 and 3 
applies, without regard to the homogeneous solution to the in- 
tegral equation. For type (1) an additional analysis similar to 
that given here must be used. 
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On Plastic Buckling of a Compressed Strip 


P. Cicala 
Instituto Aerotecnico and University of Cordoba, Argentina 
February, 1950 


1949, 
529), Professor Bijlaard came to the conclusion 


N AN IMPORTANT COMPREHENSIVE PAPER (September, 

JOURNAL, p. 
that the usual incremental (or flow) theory of plasticity does not 
predict the right buckling loads of plates, whereas a perfect agree 
ment with test results is obtained by a deformation theory. In 
an interesting discussion (above issue, p. 567), Professor Drucker 
restated the inconsistency of deformation theories and pointed 
out that a more complex incremental theory than the Prandtl- 
Reuss should be used if initial imperfections of the specimens did 
not prove to be responsible for the mentioned disagreement. 
Here, it is shown that, if a small initial deflection is taken into 
account, the incremental theory can afford an explanation of ex 
perimental facts in plastic buckling of compressed plates, even 
in the case where the discrepancy appears most marked. 

Let s be the thickness and } be the width of the strip. One of 
its long sides is free; the other is simply supported. Let o be 
the compressive stress in the direction of the length L > b. 
The twist of the buckled strip produces shear stresses r and 
strains y. We suppose that the ratio r/y has the same value 
anywhere in the plate, depending only upon o. According to 
this hypothesis, if, under the stress o, the angle of twist at station 
z passes from the initial value By) to By + 8, the maximum shear 


strain at the same station attains the value s(dB/dz). Let 


yo = s(dBo/dz) (1) 


If 8) varies according to a sinusoidal halfwave, 8 also follows 
the same law; the condition of equilibrium of moments yields 
(2) 


(s?/b?) [r (y + yo) 


c= 
Supposing Poisson’s ratio =0.5, the Prandtl-Reuss equations 

give 

(G dy — dr)a (3) 


(3G de — da)r = 
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where ¢ is the longitudinal strain. Expressing the octahedra) 
shear stress as a function of the corresponding strain vields 7 


f(y? + 3e?) ' 


7? + (07/3) = ; 


The function f, a characteristic of the material, is deduce; 
In particular, if from 


shear test the tangent modulus of rigidity G, = 


from simple stress-strain measurements. 
dr/dy is oh 
tained, the derivative f’ of f with respect to its argument is 


ft’ = (re/Ye)Gt 


where 7, and y, are the ‘“‘equivalent”’ stress and strain 


To = 


y é 
V r? + (¢7/3); ye = VV y? + 3e f 
From Eggs. (2), (3), and (4), we get 


do (: 1 “) n € s? a ; bh? 
=] 7 — jf - T 7 
dy \3 o : y¥ + vo \0? G/ | aM 


The integration of this equation may be performed by step-by. 
In the here, the 
accuracy was improved by successive approximations; at each 
for the 


step procedures. calculations carried out 


step, new value of y, the quantities were computed 


according to the following cyclic order: 


do do 
(emess); ¢; 7; Fes Ves J 5 

ay dy 

If deformation theory is used, Eq. (3) must be replaced by 


jer = Yo 8 
In this case the equations reduce to a simpler form; they may 
be solved directly for a given value of y. 
The present solution is approximate, since the ratio r/y is 
supposed to take the same values all through the plate.* In 


* Similarly, in the study of plastic buckling of columns, the ratio of the 
difference between local and mean stress to the difference between local and 
mean strain may be supposed to be constant over the column length; this 
assumption leads to the same equations as the idealized scheme used by 
Shanley (May, 1947, JouRNAL). In a paper to be published by the writer, 
it is shown that the results of the rigorous solution do not differ substantially 
from those obtained by Shanley 
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order to relate the quantities that appear in the above equations 
to the deformation of the strip, we must define the point which + 
and y refer to. If these are the highest values in the plate, the 
angle of twist at the mid-section is 


Bmax. = yL/xs (9) 


The approximation principle here used might prove helpful in 
problems of plastic buckling where the rigorous solution presents 
enormous difficulties. The precision may be improved by divid- 
ing the plate into a number of regions and attributing to each of 
them, in the determination of the effect of buckling, the value of 
the stress-strain ratio calculated at a properly chosen point of the 
same region. 

The results of calculations carried out according to Eq. (7) are 
represented by solid lines in Fig. 1. The stress-strain relation 
above the proportional limit +, was assumed to be 
Ty (10) 


Gye = 2re — 


For te < Tp we have 


y = ylo/[(Gs?/b?) — aj} 


(11) 


At t = tT, the slope of the curves shows a discontinuity* as a 
consequence of the abrupt variation of G; from G to G/2. The 
same discontinuity exists on the dotted lines that were calculated 
according to the deformation theory, [Eq. (8)]. 

These calculations show that a considerable decrease in the 
maximum stress attainable (points M in Fig. 1) results from a 
small initial deformation. The difference between this max 
imum load and the initial buckling load for the perfectly flat 
strip (y = yo = 0) is relatively more relevant in the case of the 
curves calculated on the base of the simplest incremental theory, 
which apparently can explain the relatively low buckling stresses 
obtained experimentally. 

* If a stress-strain relation giving a continuous variation of G; were used 
this discontinuity would not appear; furthermore, the curves would show a 


more marked maximum, with rapidly decreasing values thereafter. 


A Note on the Subsonic Compressible Flow 
About Airfoils in a Cascade 


Henry W. Woolard 
Cornell Aeronautical Laboratory, Inc., Buffalo, N.Y. 
March 15, 1950 


INTRODUCTION 


Bir RECENT TREND TOWARD the use of the turbomachine in 
aircraft and in other installations has stimulated interest in 
methods for predicting the aerodynamic characteristics of airfoils 
ina cascade arrangement. Asa result, many papers dealing with 
the incompressible flow about airfoil cascades have been pub- 
lished! Although most of these are rather laborious to apply, 
the problem can be considered as essentially solved. On the 
other hand, the effect of compressibility seems to have been 
neglected in the literature except for reference 2, which is applic- 
able for only one blade angle (8 = 90°). While the technique of 
applying the Prandtl-Glauert correction*® 4 to linearized subsonic 
compressible flows is well known, the apparent lack of its applica- 
tion to cascade flows appears to justify a note on the subject. 


SYMBOLS 


*, ¥ = rectangular coordinates, fixed with respect to the mean 
stream direction (see Fig. 1) 
¢ = airfoil chord 


° 
d = distance between airfoils in cascade, measured along the 


stagger line 
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h = distance between airfoils incascade, measured perpendic- 
ular to the mean-stream direction 

o = cascade solidity, c/d 

a = angle of attack 

8 = blade angle (see Fig. 1) 

y = ratio of specific heats 

V = resultant velocity 

“« = perturbation velocity parallel to the x-axis 

v = perturbation velocity parallel to the y-axis 

r = angle between velocity and normal axis (see Fig. 1) 

M = Mach Number 

2 =V1—-— M? 

p = static pressure 

p = mass density 

q = dynamic pressure 

c. = section lift coefficient based on mean-stream Conditions 

Cp = pressure coefficient based on mean-stream conditions, 
(p — po)/qo 

Subscripts 

0 = mean-stream conditions 

| = entrance-stream conditions (i.e., conditions an infinite 
distance upstream of the cascade) 

2 = exit-stream conditions (i.e., conditions an infinite dis- 
stance downstream of the cascade ) 

LE = airfoil leading edge 

TE = airfoil trailing edge 

Superscripts 
‘ = compressible flow about a cascade 


Unprimed terms refer to incompressible flow about a cascade 


COMPRESSIBILITY CORRECTIONS 


The linearized partial differential equation for the perturbation 
potential @’ of a compressible flow may be transformed to La- 
place’s differential equation for the perturbation potential ¢ 
of an incompressible flow by means of the Prandtl-Glauert rela- 


tions* 4 
r= x’, y = Qy’ (1) 
o = ko’ (2) 
where 2 = Y1 — M2, k is an arbitrary constant, and the primes 


indicate compressible flow. This transformation permits the 
determination of a linearized subsonic compressible flow in 


the x’-y’ plane from a known incompressible flow in the x-y 


plane. 

It should be observed that the Mach Number used in Eqs. (1) 
is based on the mean-stream velocity vector—i.e., a velocity vec- 
tor that is the vectorial mean of the velocity vectors at plus and 
minus infinity (see Fig. lb). This selection of a reference velocity 
follows from the well-known fact in incompressible cascade theory 
that the mean-stream velocity must be used if the Kutta-Jou- 
kowsky theorem of lift is to apply. 

The perturbation velocities at affinely corresponding points 


in the two planes are:;* 4 


ku’ / io 


“= 
(3) 


vy = (k/Q)’ | 
Now let us assume that the airfoil shapes in the two planes are 


given by 


Ves g[x — mh cos (8B — a)| + mh 
Vin? =flx’ — mh’ cos (B’ — ao’)| + mh’ 
where 
—(c/2) < [x — mh cos(B — ao)] < (c/2) 
—(c/2) < [x’ — mh’ cos (B’ — ao’] < (€/2) 
and m = 0, «1, #2,...., The subscript m refers to the 


mth airfoil (see Fig. 1) and c is the airfoil chord 
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il shape in both planes.3. 4 
plane instead of « 
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ll 


dy" /dx = df/dx’ = 


It is also apparent that 
But from Eq. (3) 


The relation for the ances between airfoils as given 
It follows from Eqs. (4) and (5) that (also see references 3 and 4) 
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(9) and the geometry of Fig. 1 we obtain the 


Considering Eq. 7 
ng relations for the blade angles and solidities: 


followi 
cot (8B — ao) = (1/2) cot (B’ — av’) 


ir 

| 8B = cot~![(1/Q) cot (B’ — av’)] + av’ (10) 

Also, 

d=d'2 = _ lies ~) = d’sin(p’ — a’)V 2 + cot? (B’ — av’) 
sin (8B — a 

hence, 


, 


a 
= . (11) 


sin (B’ — ao’)V 2? + cot*(B’ — ay’) 
By virtue of Eqs. (7), the surface pressure coefficients and the 
lift coefficients are given by the following equations: 


C,'(B’, o, ao’, Mo) = (Cp/2)(B, @, a’) (12) 


cr'(B’, a’, ao’, Mo) = (c1/2)(B, @, ao’) (13) 


In accordance with Eqs. (42) and (13), the surface pressure 
coefficients and lift coefficient in compressible flow of a cascade 
of airfoils having specified values of 8’, o’, ay’, and M/) may be 
obtained by applying the factor 1/Q2 to the corresponding quanti- 
ties (ie., Cp and c;) for an incompressible cascade of the same 
airfoil profile at the same angle of attack but having values of 8 


and o given by Eqs. (10) and (11). 


The development of cascade theories by necessity depends 
upon the use of the mean stream as a reference; however, in the 
application to compressor design, the entrance-stream and exit- 
stream conditions are usually of greater practical utility. For 
this reason it is of interest to determine the aerodynamic rela- 
tions at infinite distances upstream and downstream of the cas- 
cade. 


It is known in the incompressible case that Vo is the vectorial 


mean of V; and V2. This amounts to requiring that 4; = —1 
and 2; = —v%. Applying Eqs. (7) it follows that u,’ = —w2’ and 
»,' = —w’. Hence (remembering that V)’ = Vo), it is verified 


, 


that V)’ is the vectorial mean of V,’ and V.’ in a manner similar 


to that for incompressible flow. The main difference being that, 
in the case of incompressible flow, the line AB (see Fig. 1) is 
perpendicular to the cascade normal axis NO (corresponding to 
the compressor axis), while for compressible flow it is not per- 


pendicular, This is due to the continuity requirement 


oi’ Vi’ cos i _ ps’ Ve’ cos e 


from which it is seen that AB is perpendicular to NO only when 


po; = po’—i.e., incompressible flow. 


The velocities V,’ and V2’ are given by 


(14) 


| 
yy 
+ 
+ 


VolV (1 — (u1'/Vo')]2 + (11’/ Vo’)? 


ll 


However, in order to be consistent with the assumptions used 
in linearizing the potential equation, we must continue to neglect 
squares and products of the perturbation velocities. On this 
basis and taking into account Eqs. (7), the following relations 
are obtained: Eqs. (14) become 


Vi’ = Vo' [1 + (1/2)(u1/Vo)] | 


i ae AED (15) 
V2’ = Vo’ fl — (1/2)(u Vo) § 

The angles between the vectors Vo’, V;', and Vo’ are 
ay’ — ao’ = ao’ — a’ = 0,/Vo (16) 


The turning angle (a;’ — a2’) is 


ay’ — ae’ = 2u,/ Vi (17) 
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The linearized Bernoulli equation® yields 


2 uy - 1 
M? = M?]1+=— (14 Me 
2 Vo 2 
: (18) 
2 — 1 
Mio MiH4.:5 (: +7 M2 
2 Vo 2 


In a manner similar to that used in deriving Eqs. (18), it is easily 
shown that 
po’ /py’ = 1 + (My2/2\ u,/ Vo) 
, , 9 , { ( 19) 
po /p2” = 1 — (Mo?/2)(1,/ Vo) f 


From Eqs. (15) and (19), the dynamic pressure ratios become 


0 1 u 
fot +--e~ my 
71 Q Vo 
(20) 
qo’ 1 uy 
- = 1 + (2 — M,?) 
q2 Q Vo 


Eqs. (20) are useful for converting pressure or force coefficients 
based on mean-stream conditions to coefficients based on en- 
trance- or exit-stream conditions. 
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The Surface Pressure Coefficient for 
Approximate Compressible Flow Solutions 


E. V. Laitone 

Associate Professor of Mechanical Engineering, University of 
California at Berkeley 

Received March, 1950 


SUMMARY 


It is shown that the expressions for the pressure coefficient must be ex- 
panded as a power series in applying either the Rayleigh-Janzen or Prandtl 
iteration methods. The previous investigations have used too many terms 
for the pressure coefficient in the Rayleigh-Janzen method so that the 
maximum negative value of the surface pressure coefficient was underesti 
mated. It is again pointed out, however, that an insufficient number of 
terms have generally been used in the linearized first-order Glauert-Prandtl 
solution for the surface pressure coefficient on a body of revolution. In 
this case the square of the radial perturbation velocity must be retained in 


the calculation of the surface pressure coefficient 


RAYLEIGH-JANZEN METHOD 


A‘ SHOWN IN REFERENCE | OR 2, the Rayleigh-Janzen method 
consists of an expansion of the compressible flow velocity 
potential as a power series in terms of the free-stream Mach 
Number M.., which is considered to be a small quantity. The 
first term @o corresponds to the incompressible flow solution so 
that the compressibility effects are given by 


@ = do + Modi + Matte +... (1) 


With this relation, the ratio of the square of the local velocity 
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magnitude (V) to the square of the free-stream velocity (1.) 


may be expressed as 


r\? Vo\? Vi \* V2 \? 
= + M,? + Mo‘ ee 
Ux We ho 


where (see references 1 or 2) 


(Vo/Uw)? = go,? + o,? 
(Vi/tto)? = 2Wodiz + 2ho, 1, 
2 = gi? + gi,” + 260,62, + 2o,¢2, 


(V2 ux 


sible flow. 
Cry is (1 ioc: 


By retaining the terms of order /~* we obtain the pressure coefficient that should be used with the Rayleigh-Janzen first (1/..) approxi. 


mation. 


? Vo? + Mo?Vi? M x? 
~~ = Uo* . 4 


Similarly, the pressure coefficient that should be used with the R: 


S J Vo? + Mao?2Vi? + Ma'V2? M ~? V;? \2 
Cy, = * ad 7 : = P iad 
2 ) Uo? 4 ) 


\. . wat 


The first and second Rayleigh-Janzen approximations for the 
maximum velocity about a circular cylinder are given in refer- 
ence 2, wherein the surface pressure coefficients are computed by 
Eq. (3). As shown in Fig. 1, the first approximation surface pres- 
sure coefficients computed in this manner have a magnitude less 
than that predicted by the Glauert-Prandt] (1/V1 — M..2) 
linearized theory at the higher Mach Numbers. It can be shown 
that the magnitudes are even less when the first two terms of 
Eq. (4) are used to compute the first approximation pressure 
coefficient, as in reference 3. 


However, when the first approximation velocities are sub 
stituted into Eq. (6), the correct first approximation pressure 
coefficient, then a more reasonable variation is obtained, as shown 

The correct first approximation is only slightly 
the second approximation and does not have the 
Eq. (3) 


in Fig. 1. 
less than 
anomolous decreasing divergence obtained by using 
or (4). 


As would be expected, Fig. 1 shows that the error involved in 
using the incorrect pressure coefficient equation is not so great 
for the second approximation. In this case, Eq. (3) or the first 
three terms of Eq. (4) give practically the same surface pressure 
coefficient, which is still less, however, than the magnitude ob- 
tained by using Eq. (7), the correct second approximation for 
the surface pressure coefficient. 


In reference 3, the Rayleigh-Janzen first approximation for 
an elliptic cylinder has the surface pressure coefficient computed 
by the first two terms of Eq. (4). 
siderably less than those predicted by the Glauert-Prandtl 


These magnitudes are con- 
(1/V1 — M.?) linearized theory. However, if Eq. (6) is used 
to determine the correct first approximation for the surface 
pressure coefficient, then a better agreement is obtained with the 
Glauert-Prandtl theory. 
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’ M.‘ ( yz \3 
+ 2->7y11- : 
a4 Me :) ' | ‘ 


If we substitute Eq. (2) into Eq. (4) and retain only the terms not containing 17, we obtain the pressure coefficient for incompres 
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The exact expression for the pressure coefficient C, is 


pb — Pa 2 


Cy = = : 
‘ (1/2) pot? yM a? 


fei sy 
‘ . a Mao?{1—-— — e. 
\ 2 th co? "7 


Now Eq. (3) should not be used directly with Eq. (2) because j 
essentially contains higher powers of M.. In Previous develop 
ments, the method of expanding Eq. (3) had been incorrect 
retaining too many terms. The expansion of the pressure , 
efficient in terms of M.q is : 


(Vo?/tt @?) | 


iyleigh-Janzen second (1/..‘) approximation is given by 


Vo? m\7].: 
2(1-—\(mM. + 

Uaw* Uc « 

V2 - 


It is seen that in each case the additional terms implicitly con 
tained in Eqs. (3) or (4) lead to underestimating the magnitude 
of the surface pressure coefficient at the higher Mach Numbers 
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PRANDTL ITERATION METHOD increments are considered to be small everywhere, so that, in 


In thi method, the velocity potential is given by the super- general, their squares and higher orders may be neglected. Then 
n this me , —o . 

position of extremely small perturbations on a uniform free- 
The perturbation potential and resulting velocity should be written as 


the corresponding expression for the surface pressure coefficient, 


stream flow. 


where @; is the small perturbation velocity component in the Cp = —2(b2/Uw) — (by/Uaw)? + O(y' In y (10) 
»e-ct res irection, while is the component normal to i , ? , m : 
free-stream (Uc } ener ty hi Ast . oe - fi It is important to note that the series expansion in Eq. (8) 
y is < sable to either subsonic or supersonic flow, . . . : . 
Eq. (8) is applica . = = I continues with increasing powers of /7., so that Eqs. (9) and (10) 
cease to be valid at extremely large values of M.. unless the 
perturbation velocities are correspondingly reduced. 


Uz. 
but it must be noted that the series expansion is not valid for rel- 


atively large values of 7. unless the product (¢M.) is bounded. 
For plane two-dimensional flow, it is easily shown that ¢, and 
$, are both of the order y related to the thickness of the profile, 
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However, for the symmetrical flow about a body of revolution 
it has been shown that, for either incompressible, subsonic,* or 
supersonic® flow, ¢z is of the order y? log y, while ¢y, is still of the 
order y, where y is now the radius of the body of revolution. 
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Second Approximation to Supersonic Conical Flows 
(Continued from page 334) 
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